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AUTOMATIC  CURVE  FITTING  FOR  INTERACTIVE  DISPLAY 

Won  Lyang  Chung,  Ph.D. 
Department  of  Computer  Science 
University  of  Illinois  at  Urbana-Champaign,  1975 

This  thesis  presents  a  simple,  economical,  and  reliable  auto- 
matic graphical  data  compression  algorithm  for  interactively  displaying 
curves  on  the  screen  of  CRT-like  graphics  terminals.  An  adaptive  local 
scheme  based  on  interpolating  piecewise  cubic  polynomials  with  contin- 
uous slopes  is  designed  and  implemented.  The  scheme  accepts  a  large 
number  of  data  points  as  they  are  generated  one  point  by  one  by  an  or- 
dinary differential  equations  package  and  compresses  the  data  into  a 
compact  picture  representation.  The  important  features  of  this  algorithm 
can  be  found  in  its  adaptive  use  of  two  quantitative  measures  of  approxi- 
mation—a local  least-squares  error  norm  and  a  pseudo  distance  norm,  and 
its  ability  to  regenerate  aesthetically  pleasing  curves  using  the  piece- 
wise  cubic  polynomial  interpolants  with  relatively  few  knots  by  one-sided, 
local  computational  procedures.  The  advantages  of  one-sided,  local  pro- 
cedures are  the  economical  computation  based  on  local  data  and  the  fast 
search  for  knots.  The  problem  of  numerical  stability  imposed  by  one- 
sided, local  procedures  is  solved  by  the  idea  of  extended  intervals 
which  were  intuitively  developed  and  proved  to  have  interesting  mathe- 
matical and  computational  effects  on  the  stability  and  error  behavior  of 
the  algorithm.  The  scheme  can  be  used  for  compression,  transmission, 
and  display  of  graphical  data  in  both  on-line  and  off-line  environments. 
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1  INTRODUCTION 
1.1  Statement  of  Problem 

This  thesis  is  addressed  to  the  problem  of  designing  and  im- 
plementing a  simple,  economical,  accurate  and  reliable  automatic  graphical 
data  compression  scheme  for  interactively  displaying  one-dimensional 
curves  on  the  screen  of  CRT-like  graphics  terminals. 

The  primary  goals  are  the  reduction  of  storage  and  the  efficient 
processing  and  display  of  graphical  data.  However,  our  method  can  be 
used  also  for  the  direct  input  of  graphical  data  and  can  be  applied  to  a 
wide  class  of  data  fitting  problems  with  proper  modifications.  Hence  it 
may  ultimately  facilitate  natural  man-machine  interaction  through 
graphics. 

Three  important  and  readily  distinguishable  aspects  of  the 
problem  should  be  pointed  out.  First,  we  are  dealing  with  the  problem 
of  computer  graphics,  in  other  words,  the  problem  of  representing  and 
regenerating  graphical  data.  According  to  the  classification  by 
Rosenfeld  [27],  the  problem  belongs  to  the  area  of  transforming  picture 
descriptions  into  visual  images. 

Second,  we  are  interested  in  the  automatic  compression  of 
graphical  data  which  consist  of  a  large  number  of  points  produced  one 
point  at  a  time  by  an  ordinary  differential  equation  integrator.  We 
wish  to  compress  the  input  data  points  as  they  are  generated  without 
the  need  for  storing  all  the  input  data  points,  thereby  processing  the 


input  data  on  a  local  basis.  We  also  wish  to  reduce  a  large  set  of 
numerical  data  to  a  small  set  of  output  data  representing  a  picture. 
We  are  interested  in  an  automatic  scheme,  not  an  interactive  one  as  in 
the  graphics  systems  for  interactively  designing  free-form  curves  and 
surfaces  developed  by  Coons  [8],  Forrest  [14],  and  Riesenfeld  [26]. 

Finally,  we  are  interested  in  designing  a  practical  computer 
algorithm  applicable  to  a  wide  range  of  problems,  not  a  mathematical 
model  for  a  general  description  of  curves.  Mathematical  analysis  is 
used  as  a  tool  and  is  of  secondary  importance.  In  our  analysis,  we 
strive  for  consistency  rather  than  completeness. 

The  description  of  the  problem  environment  will  help  us  to 
have  a  deeper  understanding  of  the  problem  through  exposure  of  the  nature 
of  the  input/output  data  and  the  function  of  the  total  environment.  It 
will  also  serve  as  a  typical  example  of  the  applications  for  which  the 
method  is  designed.  The  environment  of  our  problem  is  a  generalized 
modeling  and  simulation  system  with  interactive  graphics  capabilities. 
The  system  is  designed  for  computer-aided  design  of  general  network 
systems  and  graphics  are  used  as  the  medium  of  man-machine  communication. 
The  system  configuration  takes  the  familiar  form  of  time-shared  graphics 
with  intelligent  satellites  [22]  and  its  detailed  description  can  be 
found  in  references  [17,  20].  Figures  1.1  and  1.2  show  block  diagrams 
describing  the  hardware  and  software  configurations  of  the  system. 

A  sample  display  generated  by  the  plot  package  which  is  re- 
sponsible for  graphical  interactions  is  shown  in  Figure  1.3.  We  can  see 
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not  only  the  curves  but  also  the  menu  containing  commands  and  a  message 
in  the  display  output  that  is  generated  by  the  plot  package  in  which 
the  automatic  graphical  data  compression  scheme  is  embedded.  The  im- 
plementation of  the  plot  package  is  described  in  detail  by  the  author 
in  reference  [25].  An  excellent  account  of  general  principles  in  the 
design  of  such  packages  is  given  by  Newman  and  Sproull  [2}2]  and  by  Foley 
and  Wallace  [12].  Therefore,  the  following  discussion  of  the  plot  pack- 
age is  rather  sketchy. 

Three  important  characteristics  of  the  plot  package  will  be 
mentioned  before  we  focus  our  attention  on  the  main  problem.  First, 
the  graphical  nature  of  the  package:  it  is  concerned  with  the  display 
of  graphical  data,  and  therefore  requires  efficient  representation  and 
processing  of  pictorial  data.  Second,  the  interactive  nature  of  the 
package:  the  need  for  real-time  response  and  for  a  language  of  dis- 
course, i.e.,  a  graphical  command  language,  are  important.  Third,  the 
applications  nature  of  the  package:  the  requirements  to  choose  a  pro- 
gramming language,  to  maintain  device- independence,  and  to  have  a 
picture  representation  appropriate  for  problem  descriptions  are  important. 

The  plot  package  is  written  in  FORTRAN  with  special  attention 
given  to  its  size  and  speed.  All  graphics  functions  are  realized  through 
subroutine  calls  to  low-level  graphics  software.  The  command  language 
is  designed  using  an  approach  similar  to  the  GULP  system  [23],  in  which 
the  syntax  is  specified  at  a  user  language  definition  stage  and  the 
semantics,  which  signify  the  translation  to  be  carried  out  once  a  par- 
ticular statement  is  recognized  (a  "statement"  is  the  coordinates  of  a 
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joystick  hit  at  a  command  in  the  menu  area  as  implemented  in  the  plot 
package),  are  expressed  simply  as  calls  to  subroutines.  This  results 
not  only  in  a  compact  representation  of  the  command  language  in  memory 
but  also  in  ease  of  programming. 

The  most  important  characteristic  of  the  plot  package  of 
primary  interest  to  us  is  the  picture  representation.  The  picture  rep- 
resentation should  be  a  compact  data  structure  which  serves  as  the  des- 
cription of  both  problem  data  and  graphical  data.  We  solved  this  problem 
using  an  intermediate  data  structure  generated  by  the  automatic  graphical 
data  compression  scheme.  This  intermediate  picture  representation 
serves  as  a  data  base  for  which  lower-level  display  files  are  generated. 
The  properties  of  this  intermediate  data  structure  are  critical  to  the 
size  and  performance  of  the  plot  package.  The  process  is  furthur  il- 
lustrated by  the  block  diagrams  of  data  and  programs  in  Figure  1.4. 

1.2  Characterization  of  Problem 

To  gain  insight  into  possible  approaches  to  the  problem,  let 
us  consider  the  nature  of  the  input  data  and  the  desirable  properties 
of  the  output  data.  The  input  data  are  produced  by  an  ordinary  dif- 
ferential equation  integration  package  DIFSUB  [16],  and  the  output  data 
are  the  intermediate  picture  representations  generated  by  our  compression 
scheme. 

The  input  data  consist  of  a  discrete  set  of  the  numerical 
values  of  the  abscissa  and  ordinates  for  single-valued  functions.  They 


are  generated  in  large  quantity  because  DIFSUB  uses  small  steps  to  main- 
tain accuracy  and  are  produced  one  point  at  a  time  as  each  integration 
step  is  performed.  The  derivatives  are  generally  unavailable  because 
the  systems  of  equations  contain  not  only  differential  equations  but 
also  purely  algebraic  ones  as  well.  However,  this  is  not  a  problem  as 
we  will  discover  later.  Finally,  the  data  are  assumed  to  be  exact.  At 
first  glance,  this  seems  to  be  rather  artificial,  but  we  can  justify 
this  on  the  grounds  that  the  data  unlike  empirical  data  are  accurate 
enough  (up  to  five  significant  digits  in  general)  not  to  require  data 
smoothing.  This  assumption  significantly  simplifies  the  knot  selection 
process  since  it  eliminates  nonlinear  optimization  for  data  fitting.  It 
is  also  numerically  sound  because  the  final  algorithm  is  stable  and  in- 
sensitive to  roundoff  error. 

It  would  be  desirable  if  the  output  picture  representation  were 
compact,  easy  to  process,  and  a  faithful  representation  of  the  input 
data.  The  exact  definition  of  a  faithful  representation  is  a  tricky 
problem  because  it  sometimes  means  a  numerically  accurate  approximation 
and  a  visually  close  one  at  other  times.  To  be  a  faithful  representa- 
tion, the  approximating  curve  should  be  aesthetically  pleasing  without 
any  noticeable  differences  from  the  input  data. 

We  also  wish  to  reduce  the  amount  of  temporary  storage  for 
the  computation,  i.e.,  buffers  to  save  input  data  which  are  necessary 
for  the  processing  to  obtain  the  solution.  This  is  consistent  with  our 
main  goal--the  reduction  of  storage. 
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Finally,  we  wish  to  do  all   the  computation—data  compression 
and  display—in  real   time.     This  requires  a  simple  and  fast  computational 
algorithm  for  data  compression  and  approximation.     This  will  also  reduce 
the  effect  of  computational  overhead  caused  by  the  data  compression  proc- 
ess to  a  negligible  degree. 

1 .3     Survey  of  Methods 

Many  methods  of  computer  representation  and  regeneration  of 
curves  have  been  developed.     While  the  origin  of  this  area  can  be  found 
in  the  traditional   theory  of  numerical  approximation,  interest  in  this 
problem  has  been  stimulated,  first  by  the  advent  of  digital  computers, 
and  second  by  the  exploitation  of  interactive  design  in  computer-aided 
design  systems  through  graphics  terminals.     Major  efforts  have  been 
directed  toward  the  development  of  computational  procedures  which  are 
simple,  economical,  accurate,  and  flexible.     In  particular,  computational 
requirements  imposed  by  on-line  environment  and  human  factors  have  not 
only  raised  many  questions  but  also  resulted  in  the  proliferation  of 
variety  of  useful   schemes.     We  will  survey  several  schemes  that  are 
relevant  to  our  problem. 

Since  the  spline  functions  developed  by  Schoenberg  [28]  draw 
the  most  attention,  we  will   discuss  several   schemes  based  on  splines  or 
piecewise    polynomials  which  are  of  interest  to  us.     A  general   discus- 
sion of  splines  in  conjunction  with  their  application  to  computer  graphics 
is  given  by  Ahlberg  [1],  Birkhoff  and  de  Boor  [5],  and  Greville  [18]. 
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The  most  significant  contribution  to  the  efficient  computation  of  splines 
based  on  the  local  basis  functions  called  B-splines  is  due  to  de  Boor  [10] 
Freeman  and  Glass  [15]  developed  a  scheme  based  upon  nonlinear  splines 
with  particular  application  to  the  digitization  of  curves  for  graphical 
input.  Their  method  computes  the  piecewise  cubic  polynomials  through 
nonlinear  minimization  and  determines  the  points  of  interpolation  or 
knots  by  a  set  of  heuristics  and  iterative  searches  on  a  global  scale. 
A  hardware/ software  system  for  on-line  graphical  display  of  curves  called 
PHASEPLOT  was  designed  by  Dertouzos  [11]  with  the  objectives  of  graphical 
data  compression  and  compact  display  file  generation.  The  scheme  im- 
plemated  in  this  system  is  characterized  by  the  hardware  generation  of  a 
special  class  of  piecewise  cubic  polynomials.  Although  the  scheme  can 
handle  two-dimensional  line  drawings  and  the  knots  and  piecewise  cubic 
polynomials  are  computed  by  a  global  procedure,  it  is  basically  a  one- 
dimensional,  local  scheme  because  two-dimensional  curves  are  treated  as 
two  sets  of  one-dimensional  segments  by  exchange  of  variables  and  the 
knots  are  determined  by  local  maxima  and  minima.   It  is  implemented  as 
a  global  scheme  because  each  piecewise  cubic  segment  is  fitted  by  back- 
ward matching.  A  minimum  number  of  knots  is  found  by  an  exhaustive 
search  through  the  piecewise  cubic  polynominal  space.  Another  piece- 
wise  cubic  polynominal  interpolation  scheme  by  Akima  [3]  deserves  our 
attention.  Akima' s  scheme  is  based  on  a  local  procedure—each  cubic 
segment  is  computed  from  local  data--and  is  attractive  because  of  its 
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computational  simplicity  and  its  ability  to  generate  smooth-looking 
curves.  It  is  useful  for  the  drawing  or  display  of  curves  when  the  knots 
are  given.  It  is  not  powerful  enough  to  compress  data  like  the  schemes 
developed  by  Dertouzos  and  Freeman  and  Glass. 

Another  class  of  schemes  that  is  \/ery   popular  and  deserves  our 
attention  is  based  on  Bezier  curves  and  surfaces  [14].  BSzier  curves 
are  based  on  Bernstein  polynomials  and  were  motivated  by  the  need  to 
find  a  mathematical  representation  for  free-form  curves  and  surfaces  and 
to  develop  an  efficient  computational  scheme  for  interactive  geometric 
design.  A  curve  segment,  in  Be*zier's  method,  is  defined  by  a  polygon. 
This  differs  from  the  class  of  methods  given  above  in  which  a  curve 
segment  is  defined  by  knots.  Using  Bezier 's  method,  one  can  find  a  best 
approximation  to  an  input  curve  by  interactively  modifying  the  defining 
polygon.  Therefore,  it  is  best  suited  for  applications  of  an  interactive 
nature.  A  discussion  of  the  mathematical  properties  of  Bezier  curves  can 
be  found  in  Forrest  [14]. 

Recently,  Riesenfeld  [26]  developed  a  new  more  general  design 
scheme  based  on  B-spline  curves  and  surfaces  by  combining  B-spline  theory 
with  algorithms  for  computing  B-splines  and  B§zier  concept  of  using  a 
polygon  to  describe  a  curve.  The  scheme  is  therefore  a  spline  general- 
ization of  the  Blzier  method  that  uses  B-splines  instead  of  Bernstein 
polynomials.  The  local  property  of  B-splines  makes  the  scheme  very 
interesting  and  useful  and  the  scheme  has  several  advantages  over  the 
Bezier  scheme  in  terms  of  computational  efficiency  and  interactive  de- 
sign applications. 
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The  last  class  of  methods  that  we  want  to  examine  is  that  based 
on  homogeneous  coordinate  geometry.  These  methods  describe  smooth 
curves  using  parametric  rational  polynominals  [2,  13]  or  conic  arcs  [4]. 
In  general,  each  curve  segment  is  represented  by  a  certain  coefficient 
matrix  which  determines  the  shape  of  the  curve.  These  matrices  are 
computed  by  solving  linear  systems,  and  different  continuity  conditions, 
minimization  criteria,  or  parametizations  lead  to  different  matrices, 
affecting  the  behavior  of  the  curve.  The  earliest  work  in  this  direc- 
tion was  done  by  Coon  [8]  in  the  United  States  and  by  Forrest  [13]  in 
the  United  Kingdom.  Ahuja's  scheme  [2]  is  based  on  cubic-spline-like 
rational  polynomials  and  controls  the  curve  shape  through  changes  in 
the  coefficient  matrices.  A  scheme  by  Cohen  and  Lee  [7]  is  based  on  a 
particular  class  of  curves  suitable  for  a  hardware  curve  generator. 
While  a  curve  is  described  by  a  relatively  large  matrix,  these  schemes 
are  convenient  for  the  manipulation  of  geometric  forms  and  especially 
perspective  images.  Recently,  Albano  [4]  was  motivated  by  the  recon- 
struction problem  of  border  lines  for  pattern  recognition  to  develop  a 
digitization  scheme  based  on  conic  arcs  and  line  segments.  The  conies 
of  best  fit  are  determined  by  least-squares  error  minimization  and  the 
resulting  linear  system  of  equations  is  solved  to  get  the  coefficients 
of  the  conies.  The  knots  are  determined  by  global  search  and  the  conic 
arcs  and  line  segments  are  adaptively  fitted  through  these  knots.  Each 
conic  is  specified  by  six  coefficients  in  addition  to  the  coordinates  of 
the  two  end  points  of  the  segment. 
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From  the  above  examination  of  various  methods  for  the  computer 
description  and  generation  of  curves,  we  can  make  some  valuable  observa- 
tions on  the  characterization  of  the  schemes.  Different  schemes  are 
based  on  different  mathematical  models  or  basis  functions—piecewise 
cubic  polynomials  or  splines  [1,  3,  5,  11,  15],  Bernstein  polynomials 
[14],  parametric  rational  polynomials  [2,  7,  8,  13],  or  conic  arcs  [4]. 
All  the  schemes  are  based  on  global  procedures  except  Akima's  [3].  Two 
schemes  [7,  11]  are  designed  for  hardware  generation  of  curve  segments. 
Some  schemes  are  useful  for  graphical  data  compression  of  one  kind  or 
another  [4,  11,  14,  15,  26].  All  but  two  schemes  [14,  26]  use  inter- 
polation through  knots.  Although  computational  efficiency  and  storage 
requirements  vary  from  one  scheme  to  another,  all  but  one  [3]  require 
either  nonlinear  minimization  or  the  solution  of  large  linear  systems 
for  computing  the  approximation  and  all  the  data  compression  schemes  use 
time-consuming  iterative  searches  through  the  solution  space  to  find  the 
optimal  set  of  knots  and  curve  segments.  A  simple-minded  conclusion  can 
be  drawn  from  these  observations  and  the  characterization  of  our  problem, 
it  is  namely  that  a  simple  adaptive  interpolation  scheme  using  a  small 
number  of  knots  based  on  a  local  procedure  would  be  the  most  desirable 
solution  to  the  automatic  graphical  data  compression  problem. 

Our  automatic  graphical  data  compression  scheme  is  based  on 
interpolating  piecewise  cubic  polynomials  with  continuous  slopes.  This 
approximation  is  chosen  because  it,  like  other  piecewise  polynomials, 
possesses  many  nice  properties  for  curve  representation  [1,  5,  18,  24, 
28,  29]  and  it  has  a  compact  representation  requiring  only  three  values 
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at  each  knot--the  abscissa,  ordinate,  and  slope.  It  is  also  easier  to 
compute  cubic  interpolants  than  higher  degree  polynomials.  Finally 
the  piecewise  polynomials  in  the  class  C  are  more  flexible  than  splines 
and  other  smoother  polynomials  and  hence  can  better  approximate  wildly 
varying  curves. 

The  interpolant  is  generated  locally,  segment  by  segment,  using 
the  following  procedure.  A  knot  is  tentatively  selected,  and  the  free 
parameter  in  the  cubic,  the  slope  at  the  right  endpoint,  is  chosen  to 
minimize  the  sum  of  the  squares  of  the  "errors"  at  the  sample  points 
within  the  subinterval.  Two  definitions  of  error  are  used  adaptively. 
On  relatively  horizontal  curve  segments,  the  error  at  a  sample  point  is 
the  vertical  difference  between  the  approximated  and  actual  values, 
while  on  nearly  vertical  segments,  the  error  is  an  approximation  to 
the  conventional  shortest  distance  between  the  data  point  and  the  inter- 
polant. In  the  first  case,  the  equation  for  the  slope  is  linear  while 
in  the  second  it  is  nonlinear.  Hence  linear  and  nonlinear  schemes  are 
used  adaptively  to  generate  the  cubic  segments.  If  the  errors  are  close 
to  the  maximum  permissible  error,  the  knot  and  the  segment  are  accepted. 
Otherwise,  a  new  knot  is  chosen  and  the  process  repeated.  This  one- 
sided, local  process  eliminates  the  time-consuming  iterative  procedures 
used  by  global  schemes  and  requires  a  minimum  of  storage  since  only  the 
data  relevant  to  the  current  segment  need  be  stored.  The  use  of  two 
notions  of  error  provides  additional  flexibility  and  makes  possible  a 
high  degree  of  compression. 
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1.4  Conventions  and  Notations 

We  will  define  and  discuss  the  conventions  and  notations  used 
in  the  mathematical  analysis  in  the  following  chapter  so  that  we  may  use 
them  without  further  explanation  when  they  occur. 

Since  the  algorithm  processes  the  input  data  point  by  point 
and  generates  the  output  data  consisting  of  the  knots  and  slopes  for 
the  piecewise  cubic  interpolants,  the  indices  in  our  notations  are  of 
special  importance.  We  will  use  "i"  as  the  index  for  the  input  data 
and  "j"  as  the  index  for  the  output  data  for  the  purpose  of  clarity 
and  readability.  Set  notation  is  used  because  it  is  convenient  to 
present  input  and  output  data  as  a  whole. 

The  input  data  are  a  sequence  of  points  denoted  by  the  cartesian 
coordinate  pairs  (t-,f. ).  Define 

T(n)  =  {tQ,  t,,  ...,  tn>,  ti  <  tk  if  i  <  k; 

F(n)  e  {fQ,  fr  ...,  y,  f.  =  fd^),  i  =  0,  1,  ...,  n. 

The  output  data  are  a  set  of  ordered  triples  denoted  by  (x., 

•J 

y-»  m.),  where  x-  is  the  value  of  the  j-th  knot,  y.  is  the  value  of  the 

J     J  J  J 

ordinate  at  x.,  and  m.  is  the  value  of  the  slope.     Define 

\J  J 

X(N)     =     {xQ,  x-j,   ...,  xN>,  xQ  <  x]  <  ...  <  xN; 

Y(N)     =     (y0,  yr  ...,  yN>,  yj  =  f(x.)t  j  =  0,  1,  ...,  N; 

M(N)     =     {m0>  m^9   ...,  m^}, 
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The  interval  of  interest  will  be  denoted  by  I  =  [xQ,x..].  The 
j-th  subinterval  will  be  denoted  by  I.  =  [x.  i>x.]  for  j  =  1,  ...,  N. 
The  input  data  points  which  are  deleted  for  the  subinterval  I .  ,  will  be 
denoted  by  the  following  notations: 

Xj  =  lHiy   ^(JJ+T  ti(j)+2s  ■•>  ^(jJ+Sj'  ^(j+l)  =  XJ+1; 

yj  =  fi(j)'  fi(J)+T  fi(J)+2'  '  '  fi(j)+Sj.'  fi(j+l)  =  yj+l; 
where  i  (•)  is  a  mapping  from  {0,1,  ...,  N}  to  {0,1,  ...,n}  and  s.  is 

J 

the  number  of  data  points  in  I.+,  excluding  endpoints. 

We  also  define  the  set  of  indices  of  data  points  for  each  I.  by 

S.   =  (i(j-l)  +  1,  i(j-l)  +  2,  ...,  i(j-l)  +  s.}. 

The  piecewise  cubic  polynomials  interpolating  through  X(N), 
Y(N)  with  slopes  M(N)  are  denoted  by  Py(x)  and  the  segment  of  Py(x) 
for  the  subinterval  I.  is  denoted  by  P^x).  Thus,  the  continuity  and 

J  J 

interpolation  conditions  of  Py(x)  are  represented  as  follows: 

Pj(x)   s  iflj^ajCx)  +  mj3j(x)  +  y,-_-,Y.j(x)  +  yj6j(x)'  J"  =  ls   •"»  N' 

where 


otj(x)    =    (Xj-x)    (X-Xj   -|)/hjt 


2/         x„2 


, » 


Sj(x)    =    (X-Xj.^'fxj-X)^ 

Yj(x)  =  (Xj-x^Zfx-Xj.,)  +  hjl/hj 

«j(x)  -  (x-Xj., )2[2(Xj-x)  ♦  h.]/hj3.  and  hj  =  Xj  -  Xj.,1 
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Py(x)  =  {Pj(x)}"     ; 


yj.,  ■  ftxj.,)  ■  Pj(Xj-i).  yd  ■  «xj.) 


W- 


mj-l     =^Pj(x)lx=x 


m. 


j-l 


fcPj(x)|         .  J-l.   ...,N; 
V  y0  =  f0'  XN  =  V  yN  "  V 
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2  ANALYSIS  OF  THE  GRAPHICAL  DATA  COMPRESSION  ALGORITHM 

In  this  chapter  we  discuss  and  analyze  the  adaptive  algorithm 
in  detail.  In  Section  2.1  we  discuss  the  motivations  behind  our  algo- 
rithm and  present  a  detailed  description  of  the  scheme.  In  Section  2.2 
we  present  a  mathematical  analysis  of  the  continuous  analogue  of  our 
scheme.  The  continuous  case  is  easy  to  analyze,  is  of  theoretical 
interest,  and  serves  to  illustrate  the  main  features  of  the  scheme. 
Finally  in  Section  2.3  we  briefly  sketch  the  analysis  of  the  discrete 
algorithm  following  the  pattern  used  in  Section  2.2.  The  mathematical 
details  are  not  presented  because  the  formulae  involved  are  rather 
messy.  Instead  we  content  ourselves  with  indicating  how  the  analysis 
of  Section  2.2  can  be  generalized  to  the  discrete  case. 

2.1  Introduction:  Adaptive  Interpolation  and  Knot  Selection 

In  this  section,  we  will  give  precise  description  of  the 
automatic  data  compression  algorithm  by  discussing  the  ideas  and  com- 
putational procedures  for  computing  the  interpolating  piecewise  cubic 
polynomials  and  selecting  knots.  The  following  particular  method  of 
computing  the  slope  m.  was  chosen  because  it  turned  out  to  be  superior 
to  other  schemes  after  we  experimented  with  various  local  schemes  [25]. 

The  linear  scheme  is  based  on  the  minimization  of  the  weighted 
L^-norm  of  the  difference  between  the  original  data  and  the  cubic  inter- 
polant, 

n EjH |  =  js  wk[f(tk)  -  Pj(tk)]2  (2.D 
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The  following  weights  w.  are  used  in  the  computer  program: 
tk  -  tk_r  k  =  i(j-l)  +  1 
W|.  =^tk+1  -  tk_r  i(j-l)  +  2  <  k  £  i(j-l)  +  Sj-1       (2.2) 

Vl  -  V  k  =  1(M)  +  Sj 

With  these  weights,  the  summation  is  the  trapezoidal  rule  for  the  inte- 
gration of  the  square  error  over  I-.  Other  weights  can  be  used  as  long 
as  they  satisfy  the  requirements  for  stability  of  the  algorithm.  This 
is  discussed  in  Section  2.3.  The  free  parameter  of  the  cubic  polynomial 

P.  is  the  slope  m.  at  the  right  endpoint  of  I..  This  slope  is  chosen  to 
j  j  j 

minimize  ||E-|L-  Setting 


3 


i-  o 


V  j2 

and  solving  for  m.  yields 

mj =  kLWk[fk^(tk> "  ViWOT 

-  Vj-^WV  -  yjfij(tk)6j(tk'^s     \tW?     (2-3) 

This  scheme  is  vulnerable  to  the  problem  of  exponential  error 
growth.     This  was  observed  experimentally  and  is  evident  in  the  math- 
ematical analysis  of  the  scheme  (see  Section  2.3).     The  reason  for  this 
instability  is  the  possible  amplification  of  the  errors  in  the  slopes 
m-.     This  process  can  occur  as  follows.     Suppose  the  data  points  are 
concentrated  near  the  left  end  points  of  the  intervals.     This  can  lead 
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to  a  bad  choice  of  a  slope,  say  m,.  Given  this  slope  m,,  the  next  cubic 

j  j 

segment  may  have  to  twist  itself  badly  to  approximate  the  input  data, 
leading  to  an  even  worse  slope  m,+-,.  If  this  process  is  continued, 
the  resulting  interpolant  may  not  be  smooth  and  require  so  many  knots 
that  it  is  useless  (see  Figure  2.1).  It  was  found,  however,  that  we  can 
avoid  this  problem  by  using  additional  "future"  data  to  modify  the  slopes 
so  that  the  interpolant  can  better  accommodate  the  future  behavior  of 
the  curve.  This  can  be  done  by  introducing  an  extended  interval ,  defined 
as  follows. 
Definition:  The  extended  interval  for  I . ,  denoted  by  AI . ,  is  defined  by 

MjE(xj'ti(j)+ej]w1thejil- 

It  is  convenient  to  define  the  set  of  indices  of  data  points  contained 
in  AI,: 

ASj  E  (i(j)  +  1,  ...,  i(j)  +  e.}. 

Hence  we  choose  m.  to  minimize  the  weighted  L^-norm  of  the  errors  for 
the  subinterval  including  the  extended  interval  I-  +  AI,, 

J         J 

l|E,||=     I      w  [f(tk)  -  P.(t  )]2 
J  L       keS,+AS.  K        J  K 

J      J 

This  leads  to  the  following  formula  for  m.: 

J 

mj  *  k£S.LsWkCfk6J(tk) "  mMaA)6j(V 
-  "mVVW  -  VA'VV^s  +AS  WV]2 

J  J 

(2.4) 
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Py 


Figure  2.1  Stability  of  the  Linear  Interpolation  Scheme 
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In  practice,  a  few  additional  points  outside  I-  are  sufficient  to  make 
the  algorithm  stable.  Thus  we  can  improve  the  performance  of  the  algo- 
rithm while  maintaining  the  local  property  of  the  algorithm.  This  is 
illustraded  in  Figure  2.1.,  Here  the  dotted  curve  is  the  interpolant 
obtained  using  one  additional  point  beyond  the  knots,  and  the  solid 
curve  is  the  interpolant  obtained  using  the  same  set  of  knots  but  no 
additional  points. 

We  now  discuss  the  second  method  for  computing  the  m.,  the 
nonlinear  scheme.  The  development  of  the  nonlinear  algorithm  was 
motivated  by  the  observation  that  the  linear  algorithm  with  extended 
intervals  performed  very   well  for  relatively  flat  segments  of  curves 
but  not  very   satisfactorily  for  nearly  vertical  segments.  The  linear 
scheme  requires  too  many  knots  for  the  faithful  approximation  of  nearly 
vertical  curve  segments.  A  close  look  at  the  linear  scheme  shows  an 
interesting  fact:  for  nearly  vertical  segments  each  term  in  ||e.|£  can  be 
very   large  although  the  interpolating  polynomial  is  visually  indistin- 
guishable from  the  original  curve.  This  implies  that  the  direct  error 
between  f(t)  and  P.(t)  is  a  good  measure  of  closeness  for  nearly  vertical 
curves.  This  is  illustrated  in  Figure  2.2.  This  leads  us  to  define 
another  measure  of  closeness,  a  pseudo  distance  denoted  by  d.  ,  by  the 
following  formula: 

dk  .  {MV  -,'Afti/i  (,5) 

This  "pseudo  distance"  approximates  the  shortest  distance  between 
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(Wi-i) 


Figure  2.2  Comparison  of  Two  Error  Criteria 
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f(t.  )  and  P.  and  is  used  because  the  true  distance  is  difficult  to 
compute. 

In  the  nonlinear  scheme,  the  sum  of  squares  of  the  d.  for 

2 
the  subinterval  I.  together  with  a  correction  term,  denoted  by  w(m.-s)  , 

J  %J 

is  minimized  with  respect  to  m-  to  compute  the  value  of  m..  The  cor- 

j  j 

rection  term,  where  s  is  the  first-order  difference  approximation  to  m . , 

s  =  ^j'fi(j-i)+s  .^/^xj'ti(j-l)+s.^  and  w  is  a  P°sitive  wei9ht  with 
values  between  0  and  1,  is  included  for  the  purpose  of  preventing  us 

from  having  a  bad  value  of  m.  which  might  cause  trouble  in  the  later 

stages.     Thus  we  select  m.  such  that 

357  C?     Kl2  +  w(nvs)2]  =  0  (2.6) 


keSj 


to  obtain 


mj  =  s'kL.{fk'Pj(tk)}["23J(tk){1  +  (p^(tk))2} 

-  2{fk-Pj(tk)}Pj(tk)3j(tk)]/[2w{l   +   (Pj(tk))2}2]  (2.7) 

This  formula  is  clearly  a  nonlinear  equation  for  m.  and  therefore  must 
be  solved  for  m.  by  an  interative  technique.  It  is  solved  by  Steffensen's 
interation  method  with  the  initial  estimate  of  m-  obtained  from  the 
linear  interpolation  scheme.  Although  it  is  not  feasible  to  give  a 
mathematical  analysis  of  the  nonlinear  scheme,  observations  from  the 
experiments  give  evidence  of  satisfactory  performance  of  this  scheme. 
We  did  not  encounter  any  difficulty  even  when  we  used  only  this  scheme 
for  all  curve  segments,  achieving  reasonalbe  compression  and  smooth 
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approximations.  The  difference  between  the  vertical  distance 
|f(t.  )-P.  (tk)|  and  the  pseudo  distance  d.  is  significant  for  nearly 
vertical  segment,  but  we  can  easily  see  that  d.  -  |f(t.  )-P.(t. ) |  when 
the  segment  is  nearly  horizontal  (see  Figure  2.2). 

As  hoped,  the  nonlinear  method  was  found  to  be  superior  to 
the  linear  one  for  nearly  vertical  curve  segments.  This  result  to- 
gether with  the  previous  one  indicates  that  one  method  is  complementary 
to  the  other.  Hence  it  is  natural  to  use  both  methods  in  an  adaptive 
manner  in  order  to  obtain  the  best  possible  overall  performance,  i.e., 
to  obtain  smooth  approximation  and  high  compression.  The  linear  method 
(2.4)  is  used  for  relatively  horizontal  segments,  and  the  nonlinear 
method  (2.7)  only  for  nearly  vertical  segments.  The  extra  computational 
efforts  the  nonlinear  scheme  requires  to  calculate  the  value  of  the 
slope  m.  are  offset  by  its  ability  to  produce  a  compact  curve  description 
for  interactive  display. 

The  knots  are  chosen  to  maximize  the  subinterval  lengths  sub- 
ject to  condition  that  ||y-Py||  £  EPS  where  ||»||  is  a  certain  norm  and  EPS  is 
a  specified  maximum  permissible  error  chosen  to  ensure  the  faithful 
approximation.  Two  norms  for  ||y-Py ||  were  chosen  since  they  were  superior 
to  other  norms  we  have  tested:  the  maximum  error  norm  for  the  sub- 
interval  scaled  by  the  maximum  value  of  |y|  when  the  linear  scheme  is 
used  and  the  square  of  the  maximum  pseudo  distance  when  the  nonlinear 
scheme  is  used.  Of  course,  the  same  result  could  be  achieved  by  using 
two  tolerances  for  the  maximum  error,  one  for  the  linear  scheme  and  its 
square  root  for  the  nonlinear  scheme. 
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As  a  summary,  we  will   give  a  concise  presentation  of  the  auto- 
matic data  compression  scheme.     Let  x.  -.  be  the  previously  determined 
knot  and  t-/-_ , \+     be  the  present  input  data  point  to  be  processed.     We 
wish  to  find  x-  and  P-(x)  which  maximize  h.  =  x-  -  x-  -.  while  ||y-Py|| 

J  J  J  J  J       ' 

*  EPS  for  [Xj.rXj]. 

[A]     Determine  the  maximum  value  of  the  first-order  approxima- 
tion to  the  slope  of  the  input  data  for  the  subinterval   by 


s  =   (  max  |t|<-tk_1|)YMAX/TMAX 


where 


and 


tke[xrti(j_1)+r] 


YMAX  =       max      |f(t.) 


TMAX  =  tn. 


If  s  is  smaller  than  the  threshold,  then  go  to  [B].  Otherwise,  go  to  [C] 
[B]  Linear  scheme:  Compute  the  slope  m-  at  t«#.  , \+  using 


(2.4).  Determine 

EMAX  =        max 

V^j-l'S'U-n+r 
Go  to  [D]. 


f(tk)"Pj{tk)'/YMAX' 


[C]  Nonlinear  scheme:  Compute  the  slope  m.  at  t. /.  ,  \+ 
using  (2.7).  Determine 


EMAX  =  max  |d.  |2 
keS.   K 

J 


Go  to  [D], 
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[D]  If  EMAX  <  EPS,  save  the  present  point;  get  the  next  input 
point,  r  =  r  +  1 ;  go  to  [A]. 

If  EMAX  =  EPS,  t. /._ ,>+  becomes  the  knot  x-  and  go  to  [E]. 
If  EMAX  >  EPS,  the  previously  saved  point  is  used  as  a  knot  and 
go  to  [E]. 

[E]  Clear  buffers;  get  next  input  points  and  repeat  the  above 
steps  until 

XN  =  V 

2.2  Analysis  of  the  Continuous  Linear  Interpolation  Scheme 

The  results  of  the  tests  of  the  adaptive  algorithm  described 
in  Section  2.1  using  a  carefully  chosen  set  of  examples  convinced  us 
that  it  is  powerful  and  useful  for  automatic  curve  fitting.  However,  a 
good  computer  algorithm  should  be  analyzed  and  clearly  understood  so 
that  its  users  may  use  it  with  confidence.  Only  the  interpolation  based 
on  the  linear  method  is  analyzed  because  the  nonlinear  method  has  re- 
sisted analysis.  Instead  of  analysis  we  rely  on  empirical  evidence 
that  the  nonlinear  method  works  well  most  of  time. 

The  linear  algorithm  based  on  the  continuous  error  norm  will 
be  defined  and  analyzed  in  a  rigorous  manner  particularly  because  it 
is  more  interesting  mathematically  and  simpler  than  the  actual  computer 
algorithm  described  in  Section  2.1.  Moreover,  analysis  of  the  continuous 
algorithm  illustrates  all  the  main  characteristics  of  the  algorithm, 
the  only  major  difference  being  the  dependence  of  the  discrete  algorithm 
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on  the  input  data  point  distribution  as  we  will  see  in  Section  2.3.  We 
will  show  in  this  section  that  the  linear  piecewise  cubic  polynomial 
interpolation  scheme  is  numerically  stable  and  uniformly  convergent  as 
the  size  of  the  subintervals  approaches  zero. 


We  first  define  the  continuous  Lp-norm  of  the  error  as 


rX 


J   [f(t)-p,(t);fdt 
xj-i 


We  select  m.  such  that 

■J 


E, 


3m. 


=  0 


(2.8) 


(2.9) 


and  we  obtain 


m. 


■  L(mj_1 , yj.^Yj)  =   [{  J   *(t)oj(t)dt 

Xj-1 


m.  , 
J-l 


-  y 


j-1 


^J     a..(t)3.(t)dt 
Xj-1   J  J 


Xj    Yi(t)B,(t)dt  -  yfi  6.(t)e,(t)dt]/[        (8,(t)}*dt 

Xj-1  Xj-1  Xj-1 


(2.10) 


where  a. (t),   $-(t),  y . ( t ) ,  and  6.(t)  are  functions  previously  defined  in 
j  j  j  j 

Section  1.5. 

Let  m.  be  the  value  of  the  slope  at  the  knot  x.  which  is  com- 
puted  by  the  same  type  of  formula  as  (2.10)  with  m.  •,  replaced  by  y.  ,, 
the  true  value  of  the  first  derivative  at  the  previous  knot.  That  is, 


mj  =  L(yM»yJ-TyJ] 


(2.11) 


Then,  the  amplification  of  the  error  caused  by  using  the  incorrect 
value  m.  -,  is  represented  by 
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where 


rW=xj(yj-rmM)  <2-12> 


\a  =    I   J     a,(x)3,(x)dx/  {0.(x)rdx  =  -  0.75  (2.13) 

J       Jx        J         J  Jx  J 


\H  Aj-1 

or  for  the  method  with  extended  interval 


xi  = 


rX.+Ah.  fx.+Ah.       o 


aj(x)3j(x)dx/l  J   J  {&.{x)rdx 


Xj-1  -Aj-1 


3 

35a:: 

=  -  0.75  {1 — L _}  (2.14) 

l-4a.+T0a?+15a? 

\)  \)  \J 

with  a.  =  Ah./(x.-x.  , )  where  Ah.  is  the  length  of  the  extended  interval 
We  will  call  A-  the  stability  parameter.  We  will  drop  the  subscript  j 

*J  '  i      ....   . 

from  A.  for  continuous  case  since  the  A-  are  equal  for  all  j. 

J  J 

Let  e.  =  y_.  -  m.  and  e.  =  y'-  -  m..  Using  (2.12),  we  obtain 

vi  J  J  J  J  J 

■j  "'"J  *  Xej-1 • 

I  

Subtracting  y.  from  both  sides  and  using  the  definitions  of  e-  and  e., 
j  j     j 

we  obtain 


e.  =  -  xej.,  +  e..  (2.15) 

This  error  difference  equation  describes  the  error  propagation  process. 

Now  suppose  e«  =  0  and  |e. |  ^  E.  In  fact,  we  will  show  in  Lemma  3.1 

_      3 

that  e_.  =  0(hm  ),  where  h  „  is  the  maximum  interval  length.  It  is 
j     max        max 

easy  to  show  that  |e.|  is  bounded  by 
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lejl  "  iM^-l  E  (2-16> 

Hence  |e.|  is  bounded  if  and  only  if  |x|  <  1.  The  algorithm  is  defined 

vJ 

to  be  stable  if  | X |  <  1 . 

It  is  simple  to  see  from  (2.14)  that  the  linear  method  with 
extended  interval  is  stable  if  0  <  a.  <  °°  for  all  j.  The  behavior  of 
X  is  illustrated  in  Figure  2.3  as  a  function  of  a.  The  following  ob- 
servations can  be  made  about  the  relationship  between  the  extended 
interval  and  stability. 

1.  X    =  -  0.75  for  a  =  0  (without  extended  interval). 

2.  lim  X  =  1  and  0  <  |x|  <  1  for  a  e  [0,+°°). 
a-*» 

3.  |X|   =0  for  a  -  0.3.  The  stability  is  maximized  by 

choosing  the  extended  interval  such  that  Ah.  =  0.3h.  for 

j      j 

each  step. 

Once  the  question  about  the  stability  of  the  method  is  answered, 
we  then  have  to  attack  the  problem  of  error  behavior.   In  other  words, 
as  the  size  of  subintervals  goes  to  zero,  does  the  piecewise  cubic 
polynomial  interpolant  Py  converge  uniformly  to  the  original  function  y? 
If  so,  how  fast  is  the  convergence?  Or,  what  is  the  order  of  approxima- 
tion? We  will  answer  these  important  questions  by  developing  a  priori 
error  bounds  in  the  following  theorems. 

The  following  main  theorems  prove  not  only  the  uniform  con- 
vergence of  the  approximation,  but  also  that  the  piecewise  cubic  poly- 
nomial interpolant  Py  is  a  fourth-order  approximation  to  y  with  respect 
to  the  L  -norm  for  sufficiently  smooth  y. 
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Figure  2.3  Stability  Parameter  of  the  Continuous  Linear  Scheme 
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We  now  start  our  analysis  with  a  classical  linear  approxima- 
tion theorem  by  Peano  [9],  which  is  the  basis  of  our  derivation. 

Peano's  Theorem 


n-4-T 

Let  L(f)  be  a  linear  functional   defined  on  space  C       [a,b]  of 


the  form 


L(f)  = 


b     n 


:(i) 


n     Jk 


(k) 


[  I    a.(x)f^(x)]  dx+    I      I  b    f  ^   v, 
i=0     n  k=0  i=l   1K  K*ik} 


where  the  a.(x)  are  piecewise  continuous  over  [a,b]  and  the  x.^  e  [a,b], 
Let  L(p)  =  0  for  all   p  e  p   .     Then,   for  all   f  e  Cn+1    [a,b], 


L(f)  ■ 


f(n+1)(t)K(t)dt 


where 


K(t)   =  ^f  Lx[(x-t)!J], 
(x-t)J  =   f(x-t)n,  X  >  t 


0,  x  <  t, 


and  the  subscript  x  on  L  means  L  operates  on  the  argument  as  function 

of  x. 

(r) 
The  first  main  theorem  on  the  pointwise  error  for  yx  ', 

r  =  0,  1,  2,  3  is  for  linear  interpolation  without  extended  interval. 
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Theorem  2.1 


Let  y  e  C4[I].  Then, 


l|y(rWr,L'TI||yW|Lhir.t.-o.i.2>i 


where 


To 

7 

384' 

T    2v^3  +  27 
'l     532 

T2 

2  +  9hr 
24   ' 

2  +  3h^ 

T   -         r 

'3     4   ' 

hr 

max  min' 

hmax  =  max  lhiU 
max   Uj^N  J 

and 

hm.  =  min  Ih.l . 

In  order  to  prove  this  theorem,  we  will  prove  the  following  lemmas.  Our 
proof  is  similar  to  that  by  Hall  [19]  for  cubic  spline  interpolations. 
The  values  of  all  definite  integrals  appearing  in  the  remainder  of  this 
chapter  have  been  checked  by  FORMAC  programs. 

The  following  lemma  establishes  an  a  priori  bound  for  m-. 

Lemma  2.1 


Let  y  e  C4[I].  Then,  for  each  x.  e  X(N), 

J 

lijl  Myj-ijl  *  kll"y(4,'lLhJL 
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Proof 


We  will  use  Py  and  P.(x)  to  denote  the  interpolants  obtained 

using  in.  instead  of  m..  By  Peano's  Theorem, 
j  j 


fX. 


(4) 


where 


R(y)  =  y  -  P,(x)  =   J  K(x,x)y^;(x)dx 

J     Jx.  , 

J-l 


K(x,x)  =  •|tRx[(x-t)^]  =  ^t{(x-t)^Px(x-t^}. 


Hence, 


fX. 


y  -  Pj(x)  = 


J  (  |r7K(x,T)}y(4)(T)  dx. 


3x 


\H 


We  can  set  j  =  1  without  loss  of  generality.  It  is  simple  to  show  that 

for  (x-x)+/3! 

3 

m"  =  (h"^>  {-(h-x)2(10h2+8hx+3x2)  +  22h4} 
1   24hD 

Since  (0-x)+/3!  =  0  and  (0-x)2/2!  =  0,  it  follows  that 
V3T(x~t)+}  =  Vi^  +  ^"(h-T)^  (x) 


Differentiating,  we  find 


fx"  Px{IF(x'T)  +  }  =  ml3l(x)  +  JT{h-T)l&\M 


Th 


us,   since  3-.(h)  =  1   and  6-j(h)  =  0,  we  find  that 

2 
i       _.  fh   (h-x)         _        m 

y,   -  P^h)  =       i-yr^-  m,}  yW(x)  dx 

ii  JQ 
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It  finally  follows  that 

2 

ly^l  slly^'lLJj^^-^idt^Hy^iiy.   q.e.d. 

The  following  lemma  establishes  our  a  priori  error  bound  for 
m.. 

Lemma  2.2 


Let  y  e  C4[I].     Then,  for  each  x,  e  X(N), 

le.l   =   ly'-m.l  *  ]^  ||y(4)||    h3     . 
1   3 '        '   J     J '       16  "J        M°°   max 


Proof 


Using  (2.16)  and  X  =  -  0.75,  we  find  that 

ie.|   <   |1-(-A)Ji    E<_l,.y(4),|    h3 
|ej'   "    '   1+X       '   L  "  16  l|y       "-"max* 

The  next  lemma,  due  to  Birkhoff  and  Priver  [6],  gives  the  error 
bounds  for  piecewise  cubic  Hermite  interpolation. 


Lemma  2.3 


For  y  e  C4[I], 


|y(r)Vr)IL  ^J|y(4)ll  hi-,  r  =  o,  1,2,  3 


r  ".*   "°°  max' 
where 


1/311 
u0       384'  yl       216'  y2       12'  y3       2s 
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and 

Hy  =  y^^U)  +  yj6(x)  +  y^.^M  +  y.j<5(x) 

Lemma  2.4 


For  y  e  C4[I], 

||Hy(r)-Py(r)|l   ^J|y(4)|l  hi":,  r  -  0.  1,  2,  3 


max 

where 

1      1       3hr         3hr 
a0  "  64*  al  '  16'  a2  "  8  '  and  a3  "  4  * 

Proof 

From  the  definitions  of  Hy  and  Py, 
Hy(r)  -  Py(r)  =  e.  ,a(r)(x)  +  e.6(r)(x). 

It  is  simple  to  show  that 

l|Hy(r)-Py(r)|L  s|||a(r)(x)|  +  |S(r)(x)||L  | 


■ej 


sbrTI^(4)"-hL 


where 

max     .      _  i      . 
4-,  b,       1,  b2  »  - 

mm 
Thus, 


V  -Ti'bl-1'b2-lTZ-andb3  =  12/hmin 


ar  =  br/16,   r  =  0,   1,   2,   3.  Q.E.D. 
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Finally,  it  is  easy  to  prove  Theorem  2.1  by  the  triangular 
inequality, 

l|y(r)Vr)IL  '  l|y(r,Vr)JL  +  llHy(r,-Py(r,IL. 

That  is, 

Tr  =  yr  +  ap,  r  =  0,   1,  2,  3. 

In  the  following  paragraphs,  we  will  show  that  the  a  priori 
error  bounds  for  the  linear  interpolation  with  extended  interval  is  an 
interesting  function  of  the  size  of  the  extended  interval.  By  showing 
the  graph  of  the  error  bound  constants  which  are  numerically  evaluated, 
we  can  readily  observe  the  remarkable  effect  of  the  extended  interval  on 
the  error  behavior.  In  practice,  the  proper  choice  of  the  extended 
interval  at  each  step  not  only  improves  the  compression  capability  of 
the  algorithm  by  reducing  the  error  for  the  given  subinterval,  but  also 
provides  us  with  a  useful  algorithm  which  is  both  stable  and  insensitive 
to  the  input  data-point  distribution. 

We  state  the  following  theorem  about  the  continuous  linear 
interpolation  with  extended  interval  which  summarizes  the  above  dis- 
cussion. 

Theorem  2.2 

Let  y  e  C4[I].  Then, 

l|y(r)-Py(r)ll„  ^  E,J|y(4»||„h^,  r-  0.1.2.3 


39 


where 

EQ  =  Me(a)/4  +  1/384,  E]  =  Me(a)  +  v^/216, 

E2  =  [1  +  6hrMe(a)]/12,  E3  =  [1  +  12h^le(a)]/2, 

and  M  (a)  is  given  in  Lemma  2.6  as  a  function  of  a. 

For  the  purpose  of  illustrating  the  magnitudes  of  the  error 
bound  constants  the  values  of  Me(a)  are  plotted  as  a  function  of  a  in 
Figure  2.4.  This  graph  shows  that  M(a)  and  therefore  the  size  of  error 
fluctuates  as  the  size  of  extended  interval  is  changed.  The  minimum  of 
M  (a)  occurs  at  a  -  0.3,  which  coincides  with  the  value  of  a  which 
minimizes  the  stability  parameter  X.  Therefore,  we  can  control  the  size 
of  the  extended  interval  to  obtain  the  best  performance  of  the  algorithm. 

For  comparison,  an  arrow  is  drawn  around  0.043  on  the  vertical 
axis  in  Figure  2.4  to  show  the  error  bound  constant  for  cubic  spline 
interpolation  which  was  obtained  by  the  same  procedures  as  used  in  this 
section  [19].  Even  though  this  cannot  be  used  as  a  sole  measure  of 
goodness  of  the  interpolation  scheme,  it  is  mentioned  just  as  one  il- 
lustration of  how  good  our  scheme  is. 

Since  the  proof  of  Theorem  2.2  is  similar  to  that  for  Theorem 
2.1,  we  will  just  state  the  following  two  lemmas  about  the  error  bounds 
for  m.  and  m..  We  assume  that  the  extended  intervals  are  a  constant 
function  of  the  intervals  I.. 

J 
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1.11 


1.48 


1.85 


Figure  2.4  Error  Bound  Constant  Me(a) 
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Lemma  2.5 


where 


Let  y  e  C4[I].     Then,  for  each  x-  e  X(N), 
\e.\-   l^-mjl   ^e(*)\\y{A)\\ymx 


fl+a  (1-t)|      _ 
Ke(a)  =  J0     1—2! ml'dT> 

m^   =  {(l+a-T)J[105a(l+a)2-21(l+a-T)+(3a2+4a+l)+  7(l+a-x)2(3a+2) 
-  3(l+a-x)3]  +  2(l-T)^(l+a)5(60a2-55a+ll)} 
/[24(l+a)5(15a2-5a+l)]. 


and 


Due  to  the  complexity  of  the  integral,  we  cannot  get  the  error 
bound  constant  K  (a)   in  a  neat  mathematical   form.     However,   the  values 
of  K  (a)  were  numerically  computed  as  a  function  of  a.     The  graph  of 
K  (a)  is  plotted  in  Figure  2.5. 

Lemma  2.6 

Let  y  e  C4[I].  Then,  for  each  x.  e  X(N), 

vi 


lejl  ■  |yj-"jl  *  m6(«)  Hy(4>IL»4x 
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Figure  2.5  Error  Bound  Constant  K  (a) 
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where 


Me(a)  =  Ke(a)/|l+X|. 
2.3  Analysis  of  the  Discrete  Linear  Interpolation  Scheme 

We  will  discuss  the  stability  and  error  problems  of  the  dis- 
crete linear  interpolation  scheme  in  a  rather  informal  manner  with 
emphasis  on  the  practical  aspects. 

The  linear  method  based  on  the  discrete  error  norm  can  be 
analyzed  in  the  same  way  as  the  linear  method  based  on  the  continuous 
error  norm.  The  stability  parameter  is  defined  in  the  same  manner 
except  the  integrals  are  replaced  by  summations. 

Ai  =  I      ViftJe.-ltJ/         I      w.{3,(tk)}2  (2.17) 

J        keS.+AS,   K  J      K     J      K     keS.+AS.   K     J      K 

This  formula  cannot  be  reduced  to  a  neat  form  like  (2.14).  However, 
we  can  easily  observe  that  A.  is  a  function  of  not  only  the  extended 
interval  size  but  also  the  distribution  of  data  points.  At  first 
thought,  this  sounds  like  a  rather  severe  limitation  on  the  actual 
algorithm  because  this  implies  that  the  algorithm  when  implemented  as 
a  computer  program  cannot  handle  a  variety  of  data  points.   If  this 
were  true,  the  algorithm  would  not  be  useful  since  the  input  data  is  not 
tailored  to  the  algorithm. 

For  the  discrete  linear  interpolation  scheme,  the  value  of 
the  stability  parameter  varies  for  each  subinterval  according  to  the 
data  point  distribution  within  the  subinterval.  Therefore,  the  def- 
inition of  stability  for  the  continuous  scheme  cannot  be  used  here.  We 
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will  define  the  stability  of  the  discrete  linear  scheme  as  follows:  if 
| A. |  <  1  for  all  j,  then  the  scheme  is  said  to  be  strongly  stable;  if 
| A. |  <  1  for  most  j,  then  it  is  said  to  be  (weakly)  stable.  This  means 
that  the  discrete  linear  interpolation  scheme  is  stable  if  |\. |  ^  1  for 
only  a  few  subintervals.  As  we  will  see  later,  only  the  discrete  linear 
interpolation  scheme  with  the  extended  interval  can  be  strongly  stable. 
Therefore,  the  extended  interval  is  critical  for  the  discrete  scheme 
whereas  it  is  used  only  to  improve  the  performance  of  the  continuous 
scheme. 

The  stability  problem  for  the  discrete  linear  interpolation 
scheme  is  discussed,  first  without  extended  intervals,  and  then  with 
extended  intervals.  For  the  first  case,  we  will  discuss  the  stability 
in  terms  of  data  point  distributions.  We  will  examine  the  behavior  of 
the  stability  parameter  A.  and  the  various  data  point  distributions 
which  make  the  algorithm  stable.  For  the  second  case,  we  will  explain 
how  we  can  obtain  stability  through  control  of  the  extended  interval. 
For  both  cases,  we  will  discuss  the  case  of  equidistant  data  points  as 
a  specific  example. 

Rewriting  the  expression  of  stability  parameter  for  the 
discrete  algorithm  from  (2.17)  and  setting  t.  =  £.h.  where  h.  =  x.  -  x.  -, 


and  0  £  L    $   1 ,  we  get 


j  j  j 
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Although  the  w.  are  given  by  (2.6)  for  practical  purposes,  we  assume 
that  they  are  just  any  positive  numbers  in  the  following  discussion. 
From  this  expression  for  X.  and  the  graphs  of  the  factors 
A.  and  B.  in  the  numerator  and  denominator,  which  are  given  in  Figure 
2.6  as  functions  of  £,  the  following  observations  can  be  made  about 
the  stability  parameter. 


1.  The  A.  are  functions  of  the  data  point  distribution, 

J 

2.  If  all  the  £k  >  (<)  ^,  then  |X.|  £  (>)  1.  Furthermore, 


lim  X.  =  -  °°  and  lim  X.  =  0. 

This  is  obvious  from  the  expression  for  X.  or  Figure  2.6 
If 


IwJa.-U.)]2  *  I   w.DUt.)]2, 
keS.  K  J  K     keS.    J  K 

J  J 

then  X.  <  1.  This  can  be  easily  proved  by  Cauchy's  in- 
equality.  It  should  be  noted  that  the  observation  2  is 
a  special  case  of  this. 

Let  wk  be  equal  and  let  g(£)  =  [a(?)]2-[3(?)]2.  If  for 
each  £.  <  2   there  is  a  £.  >  j  in  tne  interval  I.  where 
9(Cl)  ^  -  g(?,-)»  then  X.  <  1.  This  can  easily  be  seen 

K  1  J 

from  the  graph  of  g(£).  It  should  be  noted  that  this  makes 

plausible  the  assertion  that  X.  <  1  for  equidistributed 

data  points. 

If  the  £k  and  wk  are  such  that  A,  ^  A2  >  . . .  >  As  ,  and 

j 
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Figure  2.6     Data  Point  Distribution  and  Stability 
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B-j   ^  A-j,  B.jB2  *  A]A2,    ...,  6^2... Bs     *  A1A2...Ag    ,  then 
X.  <  1.     This  is  the  trivial   consequence  of  the  theorem 
in  Section  3.9.22  in  Mitrinovic  [21]. 
To  be  specific,   let's  look  at  equidistant  data  points  with 
the  w^  given  by  (2.6).     The  X.  become 

\    =      I  k3(s.+l-k)3/     I  k4(s.+l-k)2 
J       keS,         J  keS,         J 

It  is  easy  to  show,  using  Bernoulli  polynomials  and  the  integer  sums, 

that 

s^(l/4-3/5+3/6-l/7)+e(s^) 

X.   =  ^-7 7^- 

J     sj(l/5-2/6+l/7)+6(sV) 

Thus,  we  get 

lim  X.   =  -  0.75. 

The  discrete  linear  scheme  with  the  extended  interval  is  of 
particular  interest  because  it  is  the  actual  algorithm  we  use  in  the 
computer  program.  The  actual  computer  algorithm  uses  only  one  extra 
point  for  the  extended  interval  because  it  was  observed  from  the  test 
that  one  extra  point  simplifies  computation  without  affecting  the  per- 
formance of  the  algorithm.  Thus,  A.  is  written  as 

k  I  Vk3OW<1+*)3         4n+,3       n 
keSi                                                 a   (1+a)     -  D-, 
^   = J = !_ 

J  I  wk^(l-?k)2+wea2(l+a)4       a3(l+a)4  +  D2 
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where 


Di  ~-  i  *£u-h? and  D2  ■  i  \4"</- 

The  following  observations  can  be  made  about  the  stability 


parameter. 


1.  X.  is  a  monotonically  increasing  function  of  a,  for  any 

data  point  distribution.  It  can  be  easily  shown  that 

%—  X.  Z   0  since  D, ,  D0,  a  >  0. 
da  '  j  \c. 

2.  lim  X.  =   1. 
a-*»  J 

3.  lim  X.   =  lim  X.  =  tt-  for  a  >  0. 

^°J    VlJ    1a 

4.  -  °°  <  X.  <  1  for  0  <  a  <  +  °°. 

J 

5.  X.  =  0  for  only  one  value  of  a  £  0. 

This  is  an  obvious  consequence  of  the  observations. 
We  have  drawn  the  curves  of  X.  as  a  function  of  a  for  a  few 
data  point  distributions  in  Figure  2.7.  The  above  properties  of  X.  can 
be  easily  seen  from  these  curves.  Therefore,  the  use  of  one  extra  point 
can  be  justified  by  the  fact  that  an  extended  interval  of  small  size  for 
each  subinterval  is  enough  to  make  the  algorithm  stable.  In  practice, 
choosing  a  proper  size  of  the  extended  interval  for  each  step  requires 
the  time-consuming  process  of  evaluating  the  value  of  X.  for  each  sub- 

interval  and  extra  point  and  therefore  of  doubtful  cost  benefit.  We  can 

-2 
easily  show  that  |X,|  =  0  for  0  <  a  <  0.3  because  |D^|  <  0.8  x  10 

where  D-,  is  defined  above.  This  upper  bound  for  |D,  |  is  obtained  from 

the  area  of  the  minimal  triangle  surrounding  the  curve  A ^  in  Figure  2.6. 
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Figure  2.7  Stability  Parameters  for  the  Discrete  Linear  Scheme 
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Figure  2.8  Stability  Parameters  for  Equidistant  Data  Points 
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Finally,  we  will  look  at  equidistant  data  points  and  make  the 
following  observations  on  the  X.: 

a4(Ha)3  -  1 

1.  lim  X.  =  -^ 2 ^r-  E  GU)- 

2.  G(0)  =  -  0.75. 

3.  lim  G(a)  =  1. 

4.  G(a)  =  0  has  a  single  root  at  a  -  0.3. 

5.  For  a  >  >  1 ,  G(a)  ■*  yfr. 

The  graphs  of  the  X.  fors.  =  2  and  s.  =  100  are  shown  in 
Figure  2.8,  for  the  equidistant  data  points. 

Our  discussion  of  the  a  priori  error  bounds  for  the  discrete 
linear  interpolation  scheme  will  be  very   informal  and  therefore  limited 
to  the  statement  of  observations  and  experience  with  test  runs.  It  is 
easy  to  show  that  the  a  priori  error  bounds  are  expressed  in  terms  of 
the  X.  whose  values  vary  for  each  subinterval  depending  on  the  data 
point  distribution. 

Solving  the  difference  equation 


e.  =  -  A.e.  ,  +  e. 
J     J  J-l    J 


yields 


e*  ■  l      n   (-A.)e.  +  n  (-x.)en 
J   i=l  k=i+l   K  n   i=l   n  u 


—  (r)  (r) 

If  |e-|  is  bounded,  then  |e-|  is  bounded  and  the  bound  for  ||yv  -Pyv  'H, 

can  be  expressed  in  terms  of  the  bound  for  |e.|  as  was  done  in  Section 

•J 

2.2.  However,  the  X.  are  also  dependent  upon  the  size  of  the  extended 
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intervals.     This  not  only  makes  the  error  bound  expressions  very  messy 

but  also  makes  it  worthless  to  obtain  worst  case  error  bounds  since 

(r)       (r)n 
||  y      -Py       H^   <  +  °°  is  meaningless.     However,  it  is  not  difficult  to 

(r)       (r) 
observe  that  the  a  priori   bounds  of  ||yv     -Pyv      H^for  the  discrete 

scheme  with  the  extended  interval  will   look  very  much  like  the  graph 

in  Figure  2.4.     The  only  difference  would  be  the  magnitude  of  the  error 

bound  constants.     We  could  observe  this  when  we  actually  calculated  the 

error  bound  constants  for  a  given  sample  data  point  distribution.     The 

most  important  fact  here  is  that  the  role  of  the  extended  interval   is 

critical   for  the  stability  and  nice  error  behavior  of  the  discrete 

linear  interpolation  scheme. 


Finally,  we  will   show  that   |e.|   is  bounded  in  the  following 

vJ 


lemma. 
Lemma  2.7 


Let  y  e  C4[I].     Then,  for  each  x,  e  X(N), 


|e,l  =  |y;-I,l  <^l|y(4)ILhi 


'j  -ji       24  l|jr       M«Mmax 
Proof 

Set  j  =  1  without  loss  of  generality.     From  Peano's  Theorem, 
we  find  as  in  Lemma  2.1  that 


y-,  -  m1 

where 


2 
h     (h-x)        _        U) 

i-jT^-  m1]y(4)(T)dx! 
0 
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/    Iwk[6l(tk)]2 

keS, 


Since 


(h-T)J 

— 21-   -  m-j   ^  0  (see  Lemma  2.8)  and, 

fh  (h-x)J     _ 


{ 


-  m^dx  =  24, 


it  follows  that 


|yi-i|l  ^lly(4)ll.h3 


Q.E.D, 


Lemma  2.8 


(h-x)2 


m-,   ^  0  for  0  ^  t  ^  h,  where  m,is  given  in  Lemma  2.7, 


Proof 


Set  x  =  uh  and  t.    =  £.  h  to  obtain 

2!  ^  r         4  2  3-25.    (1-u)3 

2  (1-u);       keS.  k  k         k       '  5k         3! 


=  hfc{- 


5k»-«k-)        3! 


2! 


X  wk?k(1^k)2 
keS.       K         K 


} 


Examine  the  term  in  brackets.  Let 

g(C.u)  -  r^  -35-  -  J—  -J]- 
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We  wish  to  show  that 


(l-u)J 
g(5,u)  <   2,   for  0  <  £,  u  <  1 


It  is  simple  to  show  that 


35  g(e.u)  *  o, 


i.e.,  g(£,u)  is  a  montonotically  increasing  function  of  £,  and  that 

(l-u)J 
g(l,u)  =  —yj — 

Then,  the  desired  result  follows.  Q.E.D. 
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3  COMPUTATIONAL  RESULTS 

In  general,  there  do  not  exist  clearly  defined  criteria  by 
which  we  can  measure  the  usefulness  of  a  computer  program.  It  is 
reasonable  to  talk  about  how  good  it  is  or  how  useful  it  is  in  terms 
of  two  things.  The  first  is  the  behavior  of  the  computer  program  as 
it  is  implemented  and  executed  for  a  particular  application  or  for  a 
particular  computer.  The  second  is  the  performance  of  the  algorithm 
itself  which  is  rather  independent  of  the  particular  implementation  or 
the  run  time  environment  of  the  program.  For  a  certain  class  of  algor- 
ithms, especially  combinatorial  algorithms,  the  first  problem  tends  to 
be  treated  as  an  artistic  problem  which  is  determined  mainly  by  a 
programmer's  skill  or  style,  and  the  second  is  the  concern  of  an  area 
called  the  analysis  of  algorithms.  For  most  numerical  algorithms,  it 
is  a  common  practice  to  discuss  them  together.  However,  we  will  discuss 
them  in  two  separate  sections  for  the  following  two  reasons.  The  first 
reason  is  for  readability.  The  second  is  that  there  are  several  impor- 
tant problems  that  receive  further  scrutiny  for  the  implementation  of 
a  good  or  useful  program. 

3. 1  Implementation  of  the  Algorithm 

Since  the  existing  version  of  the  computer  programs  was 
written  and  used  for  testing  the  data  compression  algorithm  described 
in  the  previous  chapter  and  is  far  from  being  ready  to  be  integrated  into 
the  plot  package,  our  discussion  will  concentrate  on  the  various  problems 
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that  have  been  considered  to  be  important  for  the  implementation  of  the 
algorithm  and  for  the  future  improvement  of  the  present  programs.  We 
cannot  overemphasize  the  importance  of  the  factors  that  affect  the  size 
and  execution  speed  of  the  implemented  programs  in  view  of  the  fact 
that  the  algorithm  is  to  be  incorporated  into  a  package  of  interactive 
programs.  We  will  discuss  three  implementation  problems:  the  adaptive 
error  criteria,  the  data  structures  and  storage  organization,  and  the 
knot  prediction  problem. 

The  data  compression  algorithm  essentially  consists  of  two 
processes:  computing  the  piecewise  cubic  polynomial  interpolant  for 
a  given  subinterval  and  determining  the  location  of  the  knots.  The 
choice  of  an  interpolation  scheme  and  an  error  norm  for  each  segment 
is  made  adaptively  and  the  knot  is  chosen  to  maximize  the  length  of 
the  subinterval  while  maintaining  faithful  approximation.  Two  dif- 
ferent error  norms  are  used  as  criteria  for  faithful  approximation,  i.e., 
the  maximum  error  norm  scaled  by  the  display  size  when  the  linear  scheme 
is  used  and  the  square  of  the  maximum  pseudo  distance  when  the  non- 
linear scheme  is  used.  A  simple  criterion  for  switching  from  one  scheme 
to  the  other  is  used  in  the  program:  a  curve  segment  is  defined  to  be 
nearly  vertical  if  the  maximum  value  of  the  first-order  approximation 
to  the  slope  for  the  subinterval  scaled  by  the  display  size  exceeds 
a  threshold  value  THRESH.  Several  values  of  THRESH  were  tested  and 
THRESH  =  100  proved  to  be  most  effective.  As  we  will  see,  different 
values  of  the  maximum  permissible  error  EPS  are  used  for  different  test 
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_3 
examples  to  obtain  high  compression.  EPS  =  10   seems  to  be  a  reasonable 

choice  in  most  cases.  As  mentioned  before,  we  use  one  extra  point  for 

the  extended  interval  instead  of  choosing  a.  =  0.3  because  it  simplifies 

the  computation  without  significantly  affecting  the  performance  of  the 

algorithm.  When  there  is  no  data  point  within  a  subinterval,  we  compute 

the  value  of  m.  using  a  cubic  polynomial  interpolating  through  four 

neighboring  points, 

xj-r  xj'  S-UM'  ti(j)+2* 

The  second  implementation  problem  of  importance  is  the  data 
structures  and  the  storage  organization  for  the  computer  algorithm. 
The  final  data  structure  containing  X(N),  Y(N),  and  M(N)  generated  by 
the  data  compression  algorithm  is  a  simple  two-dimensional  array.  When 
there  is  more  than  one  curve  or  function  (this  is  the  case  for  our 
application),  we  need  to  have  more  than  one  such  array.  Then,  there 
exist  two  alternatives  for  the  organization  of  the  data  which  will  also 
affect  the  coding  of  the  program  in  a  minor  way.  We  can  have  either 
the  same  set  of  knots  for  all  functions  or  separate  sets  of  knots  for 
each  function.  The  former  scheme  is  considered  to  be  superior  to  the 
latter  due  to  the  following  two  reasons.  The  first  reason  is  that  the 
second  scheme  is  more  time-consuming  computationally  and  requires  larger 
arrays  than  the  first  one  even  though  it  can  reduce  the  number  of  knots 
because  the  operation  of  keeping  track  of  different  knots  for  each  func- 
tion becomes  very   messy  for  the  one-sided  step-by-step  computation.  The 
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second  reason  is  that  the  first  scheme  is  better  suited  for  interactive 
graphical  processing  since  the  use  of  one  set  of  knots  would  result  in 
faster  display  of  curves  and  more  efficient  searching  through  the  arrays. 
In  practice,  we  maintain  separate  arrays  for  X(N),  Y(N),  and  M(N),  a 
one-dimensional  vector  for  X(N),  a  multidimensional  array  for  the  set 
of  Y(N)'s  and  a  multidimensional  array  for  the  set  of  M(N)'s.  The 
question  of  the  size  of  these  arrays  is  discussed  together  with  the 
problem  of  temporary  storage. 

We  will  discuss  the  problem  of  how  to  manage  the  buffers  which 
provide  for  the  work  space  for  the  algorithm.  This  problem  is  important 
because  extra  storage  in  addition  to  the  final  data  structure  should  be 
set  aside  to  save  temporarily  a  small  set  of  input  data  points  relevant 
to  each  computation  step.  To  be  specific,  it  is  necessary  to  maintain 
buffers  for  T(n)  and  F(n),  i.e.,  (tw.  1  x+1 ,  ...,  t.^ ^+s   ,  t./jj, 
ti(j)+]}  and  {fi(j.1)+r  ....  fi(j_1)+s  .  fi(j),  fi(j)+1>  for  the  com- 
putation  step  involving  the  subinterval  I..  The  length  of  each  buffer 
when  one  extra  point  for  the  extended  interval  is  counted  is  (s.+2), 
where  s.  is  different  for  each  subinterval.  Although  the  observed 
average  of  s-  is  about  15,  we  could  observe  from  experiments  with  var- 
ious  buffer  schemes  of  fixed  or  variable  size  that  buffers  with  lengths 
of  30  to  40  words  are  needed  to  prevent  serious  degradation  in  perform- 
ance of  the  compression  algorithm.  In  one  scheme  with  buffers  of  fixed 
length  30,  we  used  a  set  of  heuristics  to  delete  appropriate  data  points 
from  the  buffer  whenever  it  was  full  until  we  obtained  the  desired  knot. 
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Unless  very  good  heuristics  are  designed,  there  seems  to  be  a  danger  of 
performance  degradation  in  the  form  of  unexpected  extra  ripples  and  bumps 
on  the  approximating  curves.  Therefore,  we  recommend  a  buffer  scheme 
which  does  not  delete  or  modify  input  data  points.  We  will  now  present 
the  buffer  scheme  which  has  been  implemented  in  our  programs  with 
reasonably  good  results.  In  this  scheme,  arrays  of  a  fixed  length,  say 
60  words,  are  defined  and  shared  by  the  final  data  structure  and  the 
temporary  data.  While  the  x-  and  the  y.  fill  up  the  arrays  from  the  one 
side  and  the  final  data  structure  for  X(N)  and  Y(N)  grows  in  size, 
the  size  of  the  usable  buffer  decreases  as  the  computation  proceeds.  It 
is  illustrated  by  the  diagrams  in  Figure  3.1.  Figure  3.1 -(a)  is  a 
snapshot  view  of  the  data  arrays  for  computation  of  the  j-th  knot.  In 
this  diagram,  the  left  side  of  the  arrays  is  being  filled  up  and  the 
computation  is  based  on  the  data  contained  within  the  arrows.  Figure 
3.2-(b)  shows  the  final  data  structure  after  the  compression  of  the 
input  data.  In  conclusion,  we  want  to  say  that  the  choice  of  the  buffer 
scheme  should  be  based  on  the  application  or  the  environment  and  that 
further  study  would  be  necessary  to  obtain  a  better  method. 

The  last,  yet  most  important  subject  of  our  discussion  of 
implementation  problems  is  the  knot  prediction  problem.  It  is  the 
problem  of  predicting  the  location  of  the  knots  by  using  certain  criteria 
or  heuristics  instead  of  just  processing  every   input  data  point.  Thereby, 
we  can  reduce  the  number  of  points  that  must  be  processed  for  the  de- 
termination of  a  knot,  and  consequently  speed  up  the  algorithm.  Im- 
plementation of  a  good  knot  prediction  scheme  would  be  very  important 
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particularly  in  a  real-time  environment.  However,  the  fact  that  the 
distribution  of  input  data  points  is  random  and  unpredictable  together 
with  the  requirement  that  the  prediction  scheme  be  able  to  handle  all 
kinds  of  curves  using  the  local  information  make  the  problem  nontrivial. 

The  prediction  scheme  can  make  use  of  two  vital  pieces  of 
information  available  at  any  given  point  of  computation,  namely  the 
past  history  of  knots  and  segments  and  the  behavior  of  the  error  norm 
for  the  present  segment.  In  particular,  the  size  of  the  error  norm 
for  a  given  subinterval  can  be  used  as  a  simple  guide  to  determine 
whether  the  present  point  is  close  to  the  actual  knot  or  not. 

For  the  purpose  of  the  systematic  development  of  a  knot  pre- 
diction scheme,  the  efficiency  of  the  prediction  scheme,  denoted  by 
E  ,  will  be  defined  as  follows: 

Lp   C  -  C  ' 
K    n    o 

where  C  ,  C  ,  and  C  denote  the  computational  efforts  with  a  prediction 

scheme, without  a  prediction  scheme,  and  with  a  perfect  prediction  scheme. 

Therefore,  the  larger  E  is,  the  better  the  prediction  scheme  is.  Since 

the  computational  effort  is  directly  proportional  to  the  number  of  input 

data  points  involved  in  the  given  computation,  we  will  use  the  number  of 

input  points  to  evaluate  the  value  of  E  for  a  prediction  scheme. 

We  will  describe  a  knot  prediction  scheme  which  was  designed 

and  implemented  with  a  reasonable  success.  It  is  based  on  the  following 

three  heuristics.  We  assume  that  we  want  to  predict  the  knot  x.+,. 

n-j+i  represents  the  i-th  predicted  subinterval  length. 
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1.  Initial  prediction:     predict  x.+l  such  that 

hil]  =  I  I     h'/Cj-D  +  hJ/2 
J+i         i=1     i  J 

The  next  subinterval  length  is  predicted  by  the  extrapola- 
tion of  the  past  subinterval  lengths.  We  put  more  weight 
on  the  immediately  preceeding  subinterval  length  h.  be- 
cause  it  is  likely  that  the  change  in  the  behavior  of 

curve  segments  is  smooth. 

(2) 

2.  Second  prediction:  predict  x.+^  such  that 


$  -  (E-^Viffl 


where  a,  is  the  error  of  the  approximation  obtained  using 
xL I   .  This  is  based  on  the  reasoning  that  the  error  is 
proportional  to  the  fourth  power  of  the  subinterval  length. 
That  is,  ^  =  K  ||y(4)IL(^+])4.  Assuming  K  and  ||y(4)|L 
are  constants  around  the  present  knot,  we  obtain  the  above 
formula. 

Prediction  using  bisection:  for  K  =  1,  2,  3,  ..., 
If  (EPS-aK+1)(EPS-aK)  >  0,  then 

h]+K|2)  =  [hg(aK+1-EPS)-h]+K;1)(VEPS)]/(aK+1-aK). 

Otherwise,  choose  a  point  whose  index  is  determined  by 

•( K+l )  (K) 

where  iv  '   and  iv  '   are  the  indices  of  the  points  which 


63 


were  selected  previously.  This  is  repeated  until  we  find 
the  desired  knot. 
The  results  of  the  experiments  with  this  knot  prediction  scheme 
are  summarized  in  Table  3.1.  Examples  1  to  5  are  described  in  Section 

3'2'  Table  3.1 

Efficiency  of  the  Knot  Prediction  Scheme 

Example  No.    No.  Input  Data  Points    No.  Knots      Efficiency  (%) 

1  209  10  74.7 

2  367  22  65.6 

3  398  6  84.8 

4  358  33  50.5 

5  120  14  60.2 

It  should  be  noted  that  different  knot  prediction  schemes  lead 
to  different  sets  of  knots  and  interpolants.  This  is  due  to  the  fact 
that  the  algorithm  does  not  exhaustively  search  through  the  solution 
space  to  find  the  optimal  set  of  knots,  i.e.,  the  minimum  number  of  knots 
which  preserves  the  smoothness  of  the  approximation.  It  has  been 
observed  from  experiments  with  various  knot  prediction  schemes  including 
"no"  prediction  that  the  choice  of  any  particular  knot  prediction  scheme 
does  not  significantly  affect  the  overall  performance  of  the  algorithm 
in  terms  of  its  compression  power  and  the  smoothness  of  approximation. 
This  is  a  very   nice  property  since  it  is  consistent  with  our  original 
objective  for  the  knot  prediction  scheme,  i.e.,  simplification  of  the 
computation  without  bad  side  effects. 
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A  source  listing  of  the  FORTRAN  programs  which  do  not  place 
a  limit  on  buffer  size  and  use  the  above  knot  prediction  scheme  is  given 
in  the  Appendix. 

3.2  Examples  and  Test  Results 

We  carefully  selected  a  set  of  five  examples  which  range  from 
a  highly  smooth  function  to  a  square-wave  function.  Although  the  tests 
were  carried  out  in  a  batch  mode  (not  on-line),  the  input  data  which  were 
generated  by  an  ODE  integrator  (we  used  a  subroutine  package  called 
DIFSUB  by  Gear  [16])  were  directly  fed  into  the  compression  package  and 
the  output  curves  were  drawn  using  a  CALCOMP  plotter.  We  would  be  able 
to  obtain  curves  of  the  same  quality  as  these  if  high-resolution  display 
terminals  were  used. 

The  test  results  are  obtained  from  the  version  of  the  computer 
programs  with  buffers  of  large  size  and  the  knot  prediction  scheme 
described  in  Section  3.1.  Computation  time  statistics  confirmed  our 
belief  that  the  data  compression  process  is  fast  enough  to  give  us  real- 
time response.  The  test  results  are  summarized  in  Table  3.2  where  the 
entries  except  for  EPS  are  the  number  of  knots  and  the  percentage  figures 
within  the  parenthesis  are  the  compression  ratios  defined  as  follows: 

r^mr^co-;™  v.^+-:«  -  Number  of  knots 

Compression  ratio  =  n— ■— —  .  ,._,,.,.  nnin+c- 
r  Number  of  input  points 

In  the  following  pages,  the  differential  equations  and/or 
circuit  diagrams  and  the  printed  output  for  X(N),  Y(N),  and  M(N)  together 
with  the  graph  of  Py(x)  are  given  for  each  example.  These  results  illus- 
trate the  power  of  the  adaptive  compression  algorithm. 
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Example  1. 

y'   =-104{y  -  sin  x  -  tanh  104(x-tt/2)}  +  cos  x  +  104  seen  104  (x-tt/2) 
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Figure  3.2 


67 


Example     2.      (Pulse  transformer): 

y'         =  7.48(l-y)  -  22.3(0.325  cos  22. 3t  -  sin  22.3t)e 

y'         =-14.96(0.0728+y)  +  44.6(0.325  cos  44.6(t-2) 

•     aa  en*  9U  -14.96(t-2) 
iin  44.6(t-2))e  , 


-7.48t 


,  t  <  2 


-  s- 


t  >  2 


y(0)    =  o 
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Figure  3.3 

X  (  J  )  Y  (  J  ) 


M(  J) 


(   1) 

0.0 

o.o 

0. 23250F 

00 

(     2  ) 

0.76911D- 

-01 

r.89993F 

00 

0.  13372F 

02 

(     3) 

0.124900 

cc 

0.1  3237F 

01 

0.27665F 

01 

(      4) 

0.15R96D 

00 

0.1.31  891= 

01 

-0.32479E 

01 

(     5) 

0.2*299D 

rr 

-•.9  3473? 

10 

-~.28170F 

01 

<     6) 

0.314310 

00 

0.90819= 

00 

0.18764F 

01 

(     7) 

0.44189D 

00 

0.10  38  3F 

01 

-0.675815 

00 

(     8) 

".514290 

00 

1.996^5= 

t)*\ 

-0.57649P 

"0 

(     9) 

0.68864D 

00 

O.1 0048E 

01 

-0.49651F- 

■07. 

(10) 

0.1RS22O 

01 

0.968Q6F 

00 

-0. 15975F 

00 

(11) 

0.20012D 

01 

0.95Q93r 

00 

-0.241?5F 

00 

(1?) 

0.2104^0 

01 

-O.1  <t640c 

oo 

0.10498E 

02 

(13) 

0.2)4870 

01 

0.44589F- 

■"1 

-0.26246E 

01 

(14) 

0.221370 

01 

-0.1  1.341  F 

00 

0.R1361P 

00 

(15) 

0.2  25310 

Ci 

-0. 7?457E- 

■oi 

0. 13650E 

CI 

(  16) 

n. 235230 

PI 

-".77767F- 

■01 

0.31202E 

on 

(  17) 

0.?446  7r> 

0? 

-0.7'  799<=- 

■01 

0.47347E- 

-01 

(18) 

0.270720 

01 

-0.7?774F- 

■01 

0.48082F- 

-01 

(19) 

0.31 758D 

CI 

-C.72B00F- 

-"I 

0.33710E- 

-CI 

(  20) 

0.3  722  4H 

0! 

-0.  7*>flC0F- 

-01 

0.26133F- 

-01 

(21  ) 

0.44433D 

01 

-0.  72800 E- 

■01 

0.29413E- 

-01 

(22) 

0.51+^60 

0! 

-0  .  72800F- 

-01 

0.29413F- 

-01 

68 


o 
ru 


o 

01 


□ 

ID 


o 

CO 


a 
o 


0,00 


0.9M 


1.8 


2.82 


3.76 


4.70 


a 

. 

i 


Figure  3.4 


Example  3.     (Fast-rise  current  switch): 


y'       =  -5.0  x  105  y  +  V/6.   x  10'5 
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V  =  0.0,  t  <  0.1 
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Figure  3.5 
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Figure  3.6 
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Example  4.     (Transformer): 


=  0.89{-161y  -  158.8cos  lOOt  +  103.542  sin  lOOt 


123(0. 16e*38t  +  1.16e"284t)} 


y(0)  =  0 


Figure  3.7 
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Figure  3.8 


Example  5.      (Double-differentiator): 

y'       =  -lOy  +  103(l-10t)e"10t 


y(o)  =  o 
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Figure  3.9 
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Figure  3.10 


Table  3.2 


Comparison  of  No.  Knots  and  Compression  Ratios 
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Example 
No. 

No.  Data 
Points 

Linear 
Scheme 

Nonlinear 
Scheme 

Adapt i  ve 
Scheme 

EPS 

1 

209 

32(15.3%) 

24(11.5%) 

10  (4.8%) 

_3 
10  J 

2 

367 

32  (8.7%) 

37(10%) 

22  (6.0%) 

2xl0"3 

3 

398 

28  (7%) 

6  (1.5%) 

6  (1.5%) 

lO"4 

4 

358 

60(16.8%) 

26  (7.3%) 

33  (9.2%) 

io-2 

5 

120 

12(10%) 

22(18.3%) 

14(11.6%) 

lO"3 
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4  CONCLUSION 
4.1  Summary 

In  this  thesis  we  have  presented  an  adaptive  local  scheme 
based  on  interpolating  piecewise  cubic  polynomials  which  automatically 
compresses  a  large  set  of  discrete  points  into  a  compact  picture  rep- 
resentation. Since  the  large  class  of  curves  representing  single-valued 
functions  can  be  handled  by  this  algorithm,  it  is  general  enough  to  find 
a  wide  range  of  practical  applications,  including  modeling  and  simula- 
tion systems,  computer-aided  design  systems,  and  statistical  analysis 
systems  with  interactive  graphics  facilities. 

The  important  features  of  this  algorithm  can  be  found  in  the 
development  and  its  adaptive  use  of  two  error  criteria  corresponding 
to  two  types  of  visual  sensitivity  to  picture  distortion,  and  its 
ability  to  regenerate  aesthetically  pleasing  curves  using  the  piecewise 
cubic  polynomial  interpolants  with  relatively  few  knots  by  simple  one- 
sided, local  computational  procedures.  The  reason  that  the  piecewise 
cubic  polynomials  with  continuous  slopes  only  can  achieve  this  much  is 
the  fact  that  they  are  more  flexible  than  smoother  piecewise  cubics 
like  splines. 

We  have  discussed  the  mathematical  properties  of  the  linear 
interpolation  scheme  together  with  the  idea  of  extended  intervals.  We 
have  shown  that  the  intutively  developed  idea  of  extended  intervals  has 
yery   interesting  mathematical  and  computational  effects  on  the  stability 
and  error  behavior  of  the  algorithm. 
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The  following  list  summarizes  the  properties  and  applications 
of  our  algorithm. 

1.  It  is  a  powerful  data  compression  scheme. 

2.  It  is  numerically  stable  and  reliable. 

3.  It  is  simple  and  therefore  computationally  economical. 

4.  It  is  flexible  and  converges  rapidly  for  wide  variety  of 
data. 

5.  It  regenerates  aesthetically  pleasing  curves. 

6.  It  has  applications  in  the  following  areas: 

a.  Compression  and  display  of  graphical  data  for  inter- 
active graphics  systems  in  which  a  large  quantity  of 
numerical  data  are  generated  by  problem  analysis  pro- 
grams , 

b.  Direct  input  of  graphical  data  or  the  digitization  of 
analog  data  for  on-line  or  off-line  use, 

c.  Compression  and  transmission  of  graphical  data, 

d.  Pictorial  data  base  generation  for  description  and 
reconstruction  of  curves  and  graphs. 

4.2  Future  Work 

We  have  mentioned  in  Section  3.1  that  further  study  and  experi- 
ments are  necessary  for  the  improvement  of  the  computer  programs.   In 
addition  to  the  development  of  a  good  buffer  scheme  and  a  knot  prediction 
scheme  which  are  moderate  in  computation  and  storage  requirements,  an 
effort  to  refine  the  code  would  be  also  helpful  to  make  the  programs 
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run  more  efficiently.  Along  with  this,  further  experiments  should  be 
carried  out  with  some  other  examples  that  contain  several  functions, 
not  just  one. 

Extension  of  the  present  data  compression  algorithm  to  general- 
ized two-dimensional  line  drawings  would  be  interesting  from  both  math- 
ematical and  practical  points  of  view.  It  could  then  be  applied  to 
such  areas  as  picture  processing  and  pattern  recognition  for  the  recon- 
struction and  description  of  pictures. 

Another  possible  direction  of  further  work  can  be  found  in  the 
extension  or  modification  of  the  present  algorithm  so  that  it  may  be 
applicable  to  more  general  problems  of  data  smoothing  and  fitting. 
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Appendix 
FORTRAN  PROGRAMS 
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DIMENSION    LBLXU  )  ,LBLY<1) 

CPMMPfsi    TFNDtYMAXfTHPESHt FPStXt XS»Yf YSfDHtNLtNRfNK 

D*T^    XLFN/5.5/,YLFN/5.5/ 

OAT*    LBLX/'  ' /,LBLY/'  » / 

HAT  4     IPtPT/1./ 

NP=4 


M* IN    PGM    FPP     AUTOMATIC    KNOT    SELECTION    ALGORITHM. 


5 

1C 
?0 
I? 

\  5 

r       dc 
r        T  * 

40 


r     c* 

•~         VK 


YM 
TH 

FH 
P  c 
F^ 
5  F 
F<~ 

F'"» 
WP 
P»" 

RF 
*C 

=  *P 

TV 

RF 
Tfv 

Tc 

TV 

xs 

YS 
TF 

LL 

=  *K 

MK 

c\ 

T  = 
DA 

Pfi 


*Xs1  .0 

PF?H= l?^. 

RM£TQM1  ,«  AOAPTT.VF.     INTERPOLATION    W/PR EDICT ICN », 5H    />LG-,II) 

AD(5,10)    ttpxT 

pvr(2r*4) 

AD (5,70)    NSV,EPS 

«m^t| T?,F10.31 

& D ( 5 , 1 2 )    Fmt 

PMAT(  18*V4) 

T"rF(6,15)     ITFXT,FPC. 

RVAT( 1HQ,2QA4/    5X , 4HFP?= , F? 0.3/ ) 

*p(5,rM'r)     rjvi(i) 

IV    INPUT    POINTS. 

AT*     pn^NTS. 

=  1. 

AD(5,FMT,FND=^0)    X(IN),Y(IN) 

=  IN+1 

c,r  yp    30 


UN. LE. 4501 

=  T  N'-l. 

(  1  )  =  X(  1  ) 

( 1  )=Y ( 1  ) 

NO=X(TM) 
A!JTnM<VTTP     KNnT 

NT'S. 

=  i 

LL    iKMHT(  TN.\'CNT,NCN,NCP) 
?#^r  v  +  mk-3 

=100.*(l.0-NCP/NCN) 

=  100.*(\'CN-NCP )/(N0N-I  ) 


SELECTION    SUP^CUTInf. 


PTTM-r    yHp    PESLH  T. 


WPITP(6,5)     M" 

WTTP(6,15)     T"TFXTtpPS 

WF  T'P(6  ,6?)     TSMKiN^F.^NtNCNTfPAtPB 
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62 


F0PM^T(5X,,«DATA=«,I5,3X,«#KNC7=«, I2,3X,»NCP=' , 15 ,3X , »NCN= ' ♦ 

♦  IS^Xt'NCNT^tl^/SXt'ABS-EFFIC^  tF6.2f  •  %  ■  ,  •  PEL-EFFIC.=  •  , 

♦  F6.  2,  *%*///     15X, «X( J)«,10X, «Y<J) • ,10X,»M(J)  •  ,10X,«H( J) •/) 
CC    7T    1=1 VNK 

IFU.LT.NK1    H«XSfl+J  )-XS(I) 
WFI7E(6, 75)     I9XS(IltYS(IltDM(  I),h 
FCRV«T(5X,M»,I2,,)'»^(2X,F12.5)) 
TTNTTNUF 


CCNTINUE 

IF(  IFLCT.LF.*>)    PETURN 

H&RDCCPY    PLOT    OUTPUT. 


160 


170 


1 
2" 


50 

C 


00  = 
X0( 
YP( 
K=l 

DC 
HX 
IF 
OX 
1  = 
OX 
XJ 
IF 
T  F 
K  = 
XP 
XI 
X2 
YP 

GC 

IF 
K  = 
XP 
YP 

cv 

CAL 

PFT 
ENC 


XS(  NK)  /3  00. 
1  )=XS(  1) 
1  )=YS(1) 


150 
=  X* 
(HX 
=  hX 
DX  + 
=  HX 
=  XP 
(  XT 
(K. 
K*l 
<K> 
=  XS 
=  XI 
(K) 

Tr 
(K. 
K*l 
<K) 
(K) 
NTT 
I  G 
UPN 


J=2tNK 
(  J  )  -X  S  (  J-l  ) 
.LE.OD)     GC    TO 
/CO 
1 

/L 

(K)+DX 

.GE.XSUH    GO 
OF. 350)    GC     TO 


170 


TC    170 
200 


=  XI 

( JJ-XT 

-XSU-1  ) 

=  (DM(  J-1)*X1*X1*X2-DM( J)*X2*X2*X1  +YS( J-l )*X1*X1* 

(2.*X2/HX*1.)*YS(J)*X2*X2*(2.*X1/HX>1.)  )/(HX*HX) 

160 
GE.350)    GC    to    200 

=  XM  J) 
=  YM  J) 

NIJF 

R*FC (xP,YP,l,Kf LBLX, 1,LBLY»1,XLEN, YLEN) 
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SUBROUTINE    ftKNOTUNtNCNT ,NCN,NCP) 
IMPLICIT    REAL*8(A-H,P-Z) 
REAL*4    Y(450),YS (50), DM< 50) 

DIMFNSIPN    X(450)  ,XS(50) 

COMMON    TEND,YMAX,THRESH,EPS,X,XS,Y,YS,DM,NL,NR,NK 

D.4TA    IWR/O/ 


C 
r 
r 
r 
f 
r 
r 


10 


20 


r 

r 

r 

30 
40 
50 


P  IECEWISE  TM 
THE  KNOTS  AR 
AND    FIXED    AN 

NCM#CrMpiJTA 
Nrp  =  «C0MPUT,« 
mcnt=  #compijt 
NT  N=o 
NCP=0 
MC  NT  =0 
TNITI  ALTZF 
IF  (IWR.NE 
+WPITP{6,3 
NL  =  1 

DIOLD=(Y( 
0ST=C!OLD 
NR=2 

DI1=(Y(NR 
D!2  =  (0U- 
ADU=DABS 
ADI2=DABS 
IFUOIl.L 
IF(ADI2.L 
I F  (  MP  .  GT  . 
IF(DIt*DI 
AOST=DABS 
IFIDMAXK 
NP  =NR  +  i 
IF(N'R.GE. 

DI0LD=0I1 

GO    tp    i  o 


BK    POLYNOMIALS    ARE    CALCULATED   ANO    THE    LOCATION    OF 
E    DETERMINED.  EACH   KNOT     IS    SUCCESSIVELY    PREDICTED 

D    THE    SLOPE    AT    EACH    KNCT    IS    CALCULATED. 

TIONS    WITHOUT    PREDICTION. 
"""IONS    WITH    PREDICTION. 
ATION-S    FOR     IDEAL    CASE. 


AND    PREDICT    THE    2-ND    KNOT. 

.0) 

10)    NKfXS(l) 

2)-Y(l) )/(X(2)-X(  1)  ) 


M  )-Y 

DIOLD 
(DID 
(012) 


(NR))/(X(NR+l)-X<NR)) 
)/(X(NP+l)-X(NR-l) ) 


E-3)    GO    TC    20 
5)    GO    TO    20 

10  TO    50 

T.    0.)     GO    TO    50 


E.    2.1 

E.    0. 

30  )     Gl 

CL.O.L" 

(DST) 

A^TI ., 4DST1/PMINICADI1, ADST  )  .GT  .100. )    GO    TO    50 


IN  )     GO    TO    50 


INITIAL    OPEPKTICN    OF    THE    NEXT    «N^T    IS    MADE. 


NL  =NP 
CALL    PRFD 
NE=N»+1 
IF(NE.GE. 
NPI-NP 
HOl=x(NP ) 
CALCULATF    tH 


KIN) 

IN)    NE=IN 

-X(NL) 

E    PIECEWISE    CUBIC    ANC   THE    ERROR    NORM 
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r 

C 

r 


TALL    PrUBE(*Ff FMAX) 
NCP  =  NCP+(NR-NL) 
IF( TWR.NE.O) 
♦  W»I"rE(6t60)     NL,NF,  X < NR) , Y <NP ) , HP1,  FMAX 
60      FORMAT<10X,2(2X,I3l t 4( 2 X , El 2 .5 ) ) 

IFUEMAX.LF.EPS)  .AND. (NR.GE. IN)  )    GC    TO    200 
TF(EVAX.EQ.ERS)    GC    TC    200 

TF(OABS(( EMAX-EPS)/EPS).LT.    0,4)    GC    TC    100 
IF  ((NR-NL.LF.l).  AND.  (  EMAX.GE.EPS))    GO    TO    200 
Ei=FV«  X 

THE    2-MC    PREDICTION'    OF    THE    KNOT. 


r*LL    PPED2(EMAX,IN,HP2) 
NE=NR+1 

IF(NF.GT.IN)     NE  =  IN 
NP2=NP 
COMPUTE    PTFCFWT^F    CUBIC    ANC   THE    ERROR, 
CALL    PCUBE(NF,FMAX) 
F2=EMAX 

NCP=NCP+(NP-NL) 
IF (TWR.NE.O) 
♦WRTTF(6,60)    NL,NP, X(NP),Y(NR) ,HP2,EMAX 
IF((EMAX.LF.FPS)  ./!MD.<NP.GE.IN)  )    GC    TO    200 
IF(FMAX.FQ.EPS)    GO    TC    2C0 

IF(  (NP.LF.NL+1) .AND. (EMJX.GE.EPS) )     GO    TO    200 
IF<( EPS-El )*<EP?-F2) .GE.O. )    GC   TO    80 
IF< MBSCNP2-MP1)     .LT.    6)    GO    TO     100 
GO    TC    9C 
TF(OABS((EPS-E?l /FPS).LT.0.25)    GO    TO    100 

THE    3-RC    PREDTCTTON    OF    THE    KNOT. 

HP  3=HP2 

CALL    PRF03(IN,HP3,E1, E2 , HP1 ,NP1,NP2 ) 

NE=N"«-1 

IF(NR.GF.IN)    NF=IN 

NP3=NP 

CALL    PCUBE(NFtEMAX) 

NCP=NFP*(NR-NL  ) 

TF UWP.NF.O) 
♦  WP TTE(6,60)    NL,NR,X<  NR)  ,Y(NP)  ,HP3,FMAX 

IF ((FM*X.LF.EPS I .AND.iNF  .GF.IN) )    GO    TO    200 

TF((FPS-Fn*(FPS-E2)     .  GT  .    0.)     GC    TF    1100 
C  IF    WF    USED    3TSFCTI0N    At G    IN    PRED3,     SEE    IF    WE 

!F(IABS(MP1-NP2)     .LE.    5)    GO    TO    100 

▼F<NP1-NP?)     1010,1010,1020 
1010       FI=F1 


BO 

C 

c 
c 

90 

91 


CAN    CONTINUE 
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1020 


1030 
1040 


1050 
1060 

r        c 
1100 


F    EPS    IS    BETWEEN    E2 


1110 


r,       i 

12  00 


C      PR 


r 
r 

r 

C 
r 

100 

110 


NPI=NP1 
FJ=F2 
N'DJ=NP? 
GO    TO    1030 
ET=E2 
NPI=NP? 
EJ  =  El 
NPJ=MPl 

IF(EMAX-FPS)     1040,200,1050 
EI =FM*X 
NPI=NR 
GO    TO    1060 
FJ=EV*X 
NO J=NP 

IF(NP.J-NP!     .IT.     6)    GC    TO     100 
OP    to    ?00  0 
HECK    IF    WE    C£M    USE    BINARY    SEARCH,     I.E.,     I 
NPI=NP2 
EI=E2 
M  P  J  =M  P  ? 

EJ  =  F2 

IF((EPS-E2)*(EDS-FMAX)     .GE.    0.)    GO    tt    1200 

IC(EMAX.GT.     F2)    GC    TC    1110 

EI=E^Ax 

N»I=NP 

GC    TP     1060 

Ej=Ewax 

No  j  =  no 

GC    ^0     1060 
F       WF    «TIt.L    CA\'    NOT    USE    BINARY    SEARCH,    CALL    PFED3    AGAIN. 

El=c2 

E2=EMAX 

NP1=NP2 

NP?=NR 

HP  1=HP2 

HP2=HP3 

GC    Tp    9" 
EDTrTTnw    BY    RISFC^ION. 

CALL     PRFDBC(EI,FJ,NPI,NPJ,NCP) 

EMAX^EI 

FUG    IS    <FT,     SHOWING   THE 
ND=-7 ,    WOVE    FORWARD. 
N0=1,    ^OVE    BACKW/iPO. 


DIRECTION    OF    rrMPUTATION 


IF(FMAX-FP^) 

IF(NF.GE.IN) 
IMD=-1 


110,200,120 

GC    TO    200 
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GO    TO    150 
120  IND=1 

IF<NP.I.E„NLH  )    GC    TQ    200 

NR=NP-1 

C* 

C   PROCEED  ONE  QY  ONE  PT  UNTIL  THE  KNOT  IS  OBTAINEC, 
r 

150    DMSV=DM(NK+1) 

TALI     PCUBE(NE,EMAX) 
NT  P*NCP*<NR-NL) 
H=X(NR)-X(NL> 
IF< IWP. NE.O) 
♦WR  ITF(6,60)    NL,NP,X(NP),Y(NR ) ,H,EMAX 
IF(E*AX-EPS)    160,200,170 

160  IF(TND)     161,161,200 

161  IF (NR.GE.TNI     GO    TO    2C0 
NR=NP+1 

NE  =  NR«-1 

IF(NF.GT.IN)     MF=IN 

IND=-1. 

GO  TP  150 
17^    IF(TMO)  250,172,172 
172    IF(NR-NL.LE.l)  GC  TO  200 

NF=NR 

*F  =NP-1 

TND  =  1 

GO  tc  150 

c 

r   PRESENT  PHIMT  IS  S*  VFD  * S  *  KNOT. 

r 

200    NK=NK-H 

YSCMK )=Y(NR ) 
XS  (NK)=X(N<?) 
GO    to    300 

r 

C       PREVIOUS    °OINT     IS    SAVED    AS    A    KNOT. 

r 

250         NK=NK+1 

NR-NR-1 

XS(NK)=X(NP) 

YS(NK  )=Y(NR  I 

DM<MK)=DM^V 


300 


PREP/^PF    PHR    thF    NEXT    KNOT. 

TNP  =  0 
I»NR-NL 
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N<"  N=  NGN  +  ( I  *  ( I  +1 )  /2-  1  I 

NCNT=NCNH-(NR-NL) 

H=X(NR)-X(NU 

I'M  IWR.NE.O) 
+  WPITf=(6,60)     NL,NP,X(NR),Y(NR)  ,H,EMAX 

IF(NP  .GE.TN)    GO    TH    320 

DX=X  (NR+1  )-X(NR) 

IF(DX.LE.    0.0)    DX=l.E-20 

P2=(Y(M0+1)-Y(MR  ))/DX 
320       DX=X(NR)-X(NR-1) 

IF1DX.LE.     0.0)    DX=1.E-2G 

ri=(Y(NP)-Y(NR-l))/DX 

IF( IWP.NF.O) 
♦  WPTT<6,310)     MK,XS(NK)  ,DM(NK)  ,01  ,02 
310       FHPM^(3X,,<«  ,12,  »>»,3X,  »X(I  )  =  •  ,E  12.  5,  3(  2X,  E  12.  5)  // 

1         13X,»NL«  ,3X,»  NR' ,6X, » X ( NP ) • ,9X, «Y(NP P,11X, • H' , 12 X, • EM AX • ) 

IF(NR.GE.IN)    RETURN 

IF('NK.GF.50)    RETURN 

GO    tp   30 

END 
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FUBPCUTINE    PCUBF(NEtEMAX) 
IMPLICIT    PFAL*8( A-H, P-Z) 
PEAL*4    Y(450)fYS (50),DM(50) 
DIMENSION    X<450)  ,XS(50) 

COMMON    TND,YMAXfTHRFSHfEPS,X,XSfY,YStDM,NL,NR,NK 
r 

C  PIECPWISE    CUBIC    POLYNOMIALS    ARE    DETERMINED    FOR    THE    INTERVALS 

C  <X(NL)«X(NR|)     BY    CALCULATING    THE    SLOPES    DM(NK+1)    AT    X(NP). 

C  AND    THFN,    ERROR    BFTV-EEN    THE    CUBIC    AND   THE    ORIGINAL    DATA    IS 

C  DETERMINEO.         TWO    DIFFERENT    FRRCR    CRITERIA    ARE    USED    ADAPTIVELY. 


■C 


NKX=NK+1 

IFINR-NL.LE.1  )    GC    TC    50 
DM<NKX)=SLPPE<NEfDMXJ 
EMX=0. 
NPTsNR-NL-1 
EMAX=0. 
H=X(NP)-X(NL) 
C       DETERMINE    THE    w»X    VALUE    OF    Y. 

n^  21  j=nl,nr 

SJ=ABS<Y< J) ) 

20  YVAX=DMAX1(YM^X,SJ  ) 

H?=H*H 
C       DFTFPMIN^    THE    FP0C^     NOPM    FCP    THE    CUBIT. 
DO    40    J = I, MPT 
K  =  fa+J 

X*=X(NP )-X(K  ) 
XP=X(K)-X(NL  ) 
XA2=XA*XA 
XB?-XB*XR 

S  J=(DM(NK)*XA2*XR-DM(NKX)*XB2*XA+Y(NL  )  *XA  2*  (  2  .*XB/H+1.  )  ♦ 
i  Y(NP)*XB2*(2.*XA/H*1. ) )/H2 

IF(CMX.GE.TH»FSH)    GO    TC     30 
FP=0*BS($ J-Y(K) J/YMAX 
EWX=^MAX1(FP,PMX) 
GC    TO    40 
30  S  JN  =  tDM<  NK  )*XA*IXA-2.*XB)-CMINKX  )  *XB  *(  2  .  *XA-XB  )  «■ 

I  6.*( Y(NR)-Y(NL) )*XA*XB/H)/H2 

EP=SJ-Y(K) 

cR=pp*cp/(j .+SJN*SJN) 
EMX=DMAXl( FR ,EMX) 
40         CCNTINUE 

EMAX=DMAX1(FMAX,EMX) 
PFTIJRN 
C 
C       NT    FXTRA    DATA     dptNT     inj    Thf    INTERVAL;       SLOPE    IS    CALCUALTED    BY 

r.       .#N     INTERPOLATING    CUBIC     POLYNOMIAL, 
r 


50 


70 


90 


NL1=M- 

fup J=NR+ 
*1  =  X(NJP 
/3=X(NR 
PI=X(NL 
IFINL.L 
£4=X(NR 
BA=X(NL 
P5=X(NL 
YD=Y(NP 
D^MNKX) 


1 
1 
I- 

)- 
)- 

E. 
>- 
1  ) 
I) 
) 
=  ( 


1 


X(NR1) 
X(NL) 

X(N»1 ) 

1)    f,n    TC    70 

XCN'U.  ) 

-X(NL) 

-X(M»1) 

Y0*B4*B5-A1*A3*Y{NL1) ) /( A4*B4*B5)  ♦{  Y0*B4*Bl-»-Al*A4* 
Y(NLM/  (A3*B4*B1)+(YC*B5*B1-A3*A4*Y(NR1I  )/(  A1*B5*B1) 


DM(MKX  )  =  Al*Y(NL)/(-A3*Bl)+<A3+Al>*Y<NR)/(A3*Al)*-A3*Y<NRl)/(Al*Bn 

RETUP^ 
FND 


91 


FUNCTIrN  SLroF(NE,DMX) 
IMPLICIT  RFAL*8(A-H,P-Z1 
PEAL*'*  YC45H) fYS (50), DM<50) 

DIMENSION  XU50),XS(50) 

rnMMCN  TEND, YMAX, THRESH, EPS, X,XS,Y, YS,DM,NL ,NP ,NK 
r 

C   THE  SLOPE  AT  X(NR)  IS  CALCULATED  FOP  THE  INTERVAL  (X (NL )  ,X ( NR )) , 
C   USING  THE  V«LUES  XS  <NK )  ,  Y(  MK  )  ,  AND  DM(NK). 
C   NK  =  pnir-Tcp  JO  THE  PREVIOUS  KNOT. 


r 


CALL    UNDERZCGFF'  ) 
XJ2=X(NR)    . 
XJl=X(NL) 
H=XJ2-XJ1 
XA=XJ2-X(MLI 
XB=X(NL)-XJ1 
BXI=-XB*X3*XA 
YPrLD=RXl*Y(NL) 
APrLC=XA*XA*XB*BXI 
BRCLD=RXI*PXT 
OBCLD=XA*XA*( 2.*XP*H)*BXI 
DPrLD=XR*XP*(  2.*XA-fH)*BXI 
YBI=0. 
ABI=C. 
PBI=0. 
CBT=0. 
DBT=C. 
NI*NL+1 
DMX=0. 
CALCULATE    th^    SLPPF    AT    X(NP),     USING    THE    LINEAR     FORMULA. 
CP    10    J=NT,*F 

DX»XU)-X( J-l) 
X4=XJ2-X( J) 
XB  =  X<  JJ-XJl 
RX!=-XP*XB*XA 
YB=PXI*Y(J ) 

AR=XA*XA*XP*BXI 

ps=BXl*nxi 

CB=XA*XA*(2.*XB*H)*3XI 

CR*XB*XB*I2.*XA*H|*BXI 

YBI=YBT  ♦(  Y3+YRPLD)*DX 

API=ABI ♦( AB+ABCLD)*DX 

RRl*B£I+(BB+BBOLD)*DX 

CBTsCBT+(CR*CBPLDt*DX 

DQT=DBH-(DB  +  DBCLD)*DX 

YPCLO=YB 

A«rLD=AB 

BBCLD=BB 
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C9CLD=CB 

DBPl  n  =  OR 

IFU.LE.NR  I    0MX=DM*X1<DMX,DABSUY(  J)-Y(  J-l)  )/DX)) 
10  r.CNTINUe 

SLC-PF=(  YBI*H*H-DK«NK)*ABI-YS(NK)*CBI/H-YCNR)*DBI/H)/BBI 

YC>=YV^X 

IF(YB.EQ.0.)     YB=1. 

HWX=CVX*TFMD/YR 

IFCDVX.LT.THRFSH)    RETURN 
C       NONLINEAR    IN  TEB&0LATI0N. 

S=(Y(NP  )-Y(W-l)  )/(X(NR|-XfNP.-l)  I 
30  SLPrLO-=SL0PE 

FX=F(XJ]  ,XJ2,S,SLP0LD) 

FFX  =  F{  XJ1  ,XJ2»S,  FX) 

PNX=PFX-?  .*FX*SLPCLD 

TF(FNX.EO.0.»    GC    TO    50 

SLCPF=5LPCLD-(FX-SLPCLC)**2/ENX 

IP(rABS(<SLr'PF-SLPOLD)/SLPCLD).LT.l  .  F-5  )    GC    TF    30 

50  cp.TUPN 

FND 


FUMCTTFN    F« XJl,XJ2fSfXK) 

IMPLICT-    P=AL*8(A-HtP-Z) 

FF«5L*4    YC+50) ,YS (50) ,CM(50) 

CTWFNSTHN    X< 450) ,XS( 50) 

CGMVCN    ^FMDtYMAXtTHPFCH.FPSTXtXStYtYS.PMtNLtNP  fNK 

FAT«    WA./ 
CMTULSTF    THF    VALUE    CF       X  =  F(X)     FQR    STEFFENSEN'S    INTERATION. 


10 


h= 

H2 
H3 
I? 
13 

or 
x 

X 

p. 
3 

P 

p 
p 


1 


F  = 
RE 
FN 


XJ2-XJ1 
=  H*H 
=H2*H 
=  "ML  +  1 
=  NP-1 
LM  =  0. 

10  1  =  1 
A=XJ2-X 
p=X(I )- 
XT=-XB* 
XT1=-XP 

i=x**rx 

Y(M 
T?=PJ l* 

(  W*f 

CMINUF 

TURN 
C 


2,13 

(I  ) 

XJ1 

X3*X 

*C2. 

H-XA 

*)*X 

K)*X 

DU 

(Y<  T 
l.+o 


A/H2 

*XA-XB) /H2 

*XA*(DM<NK)*XB  +  Y(NU*(?.*XB/H«-l.  I  )/H2  + 

B*XB*(2.*XA+H)/H3 

A*< XA-2.*XR)/H2+XK*RXI1+6.*(Y(NP)-Y<NL) )*XA*XB/H3 

)-PI)*(-BXl*(l.+PI2)-(Y( n-PI)*PIl*BXll)/ 
I?)*(l.+f>12)) 
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SUBROUTINE    PREDIUN) 
IMPLICIT    RFAL*8<A-H,P-Z) 
PEAL*4    YC450) fYS(50) ,DM<50) 
DIMENSION    X(450) ,XS( 50) 

COMMON    T END, YM AX, THRESH, EPS, X,XS,Y,YS, DM, NL,NR,NK 

r 

C       PREDTCTION    OE    THE   NEXT    KNOT    IS    MADE. 

r 

HS  =  0. 
IS  =  NK-2 

IF(TS.LE.O)    GO   TO   15 
DP    10   I "if IS 
10         HS=HS*XS (I+1)-XS< I) 
C       HS    IS    THE    SPACING    OF    THE    NEXT    INTERVAL     PREDICTED. 
15  HS=(HS/(NK-l)+XS(NK)-XS(NK-l) )*0.5 

C       INDEX    OF   THE    PREDICTED    KNOT     IS    MADE. 
IS  =  NLH 
NR  =  NL 

DO    20    T=IS,IN 
IF( X< I)-X(NL).GE.HS)     GO    TO    30 
20         CONTINUE 
30  NO=l 

IF(  <NR.LF.NL+1).AND.  (NP.LT.IN))     NR=NL+? 
IF(NP.GT.IN)    NR=IN 
RFnjPN 
END 

SUB"OUTTNE    P«ED2  (ERROR , I N, HP2 ) 
IMPLICIT    RFA|_*8(/!-H,P-Z) 
PF£L*4    Y(450),YS (50)9DM(5D 
DIMFNSION    X(450) ,XS( 50) 

COMMON    TEND, YM AX,THR ESH, EP S , X, XS, Y , YS ,DM,NL, NR , 
r 
C       THE    2-ND    PPEDKTIPM    pp    THE    KN^T    IS    MADE. 

r 

H=X(NR)-X(Nl  ) 
0  ASSUME    THE    ERROR     IS    THF    4-TH    ORDER    OF    H. 

HP2=H*(FPS/EPR°P )**0.2  5 
NROLD=NP 
IS=NL*2 
DO    10       T=I«fTN 
TF(X(I)-X(NL).GE.HP2)    GO    TO    20 
10         CONTINUE 
20  NR=T 

IF (MP.NF. NRHLD)     GO    Tp    30 
IF(HP2.GE.H)    NR=NR0LD+1 
IF(HP?.LT.H)     NR=NP0LD-1 
IF(NP.LE.NL)     NR  =  NL«-1 
30  IF(NR.GT.IN)    NR= IN 

IF((FRROP.LT.FPS*EPS).ANC.( (NR-NL) .GT .2*<NR0LD-NL ))) 
♦  NR=NL+(NR01D-NL)*2 

HP2=X(NR)-X(NL ) 
RETURN 
CND 


NK 
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SUBPPUTINF    PRED3(IN,HP3,E1,E2,HP1,NP1,NP2) 
JMP|KTT    0FAL*8(A-H,P-Z) 
PPAL*4    Y( 450) ,YS(50) tPMI 50) 
DIM-NSICN    X(450),XS(50) 

rpMMrM    TF ND, YMAX, THRESH, EPS, X,XS,Y,YS,DV,NL,NP,NK 

r 

r       THE  3-CD  PRFDITTIrN  °F  THE  NEXT  KNOT  IS  M/*9E. 

r 

~     ip    THERF    A?p    top    m^NJY    DATA    PTS    BETWEEN    PPEVICUS    PREDICTIONS, 

r       ,,SF    tHp    BISFCTION   OF    <D.*TA    pt$     fCR    THE    NEXT    PREDICTION. 

Tc(  (  EPS-F1  )*(  EPS-E2  )     .r,T.    0.)    GC    TC    1 

IF(  I*RSINDi-MD2)     .LE.    7)    GP    Tn     1 

NP  =  (Nt>l*-NP2  )  /2 

op    T    60 

1  TF{ cI.EO.p? )    GP   tp    30 

TF(( \oi-NP2)*<Pl-F2)     .GF.    0.)     GC    TC    2 
Ex  =  El 
El=c2 
F  2  =  c  T 

2  hp?^mP3 

f  LINE4P    TNTFRoiLATION    OF    PPFVICUS    2    PREDICTIONS    IS    USCD. 

md3  =  mpi+(hP3-HP!  )*(EPS-F1 ) /(E2-E1  ) 
NP?LD=NR 

IS-NL*.? 
IFIIS.GT.IN)     TS=IN 
D?    in    T=icf  tn 

H=X(T )-X(NL ) 

TF(h.LT.HD3)    00    TO     10 

TF(H-HP3.LE.HP3-X<I-1) )    GC    T0    20 

T  =  T-1 

or  -"-n   ?n 
]r         rpMTjMijF 
?0  NR=T 

IF(NR.L?.NL*1I    NP  =  NL+2 

IF ( (Fpppp .LT.EP^  *EPS  )  .ANC.<  (NP-NL ) .GT  .2*( NROLD-Nl  ) ) ) 
f  MP  =NL«-2*(NorLD-NL) 

tc (Nk .ne.ndpl D)    GC    TP    5H 

^  It ( FPS.GT.F? )     NP=NR+i 

IF(  FPC  .1  T.c?)     NR«NR-1 
5"  IF(NP#GT.INI    NR«IN 

60  HP3«X<  NP)  -X(  vIL) 

C  p  TljQM 

ENO 
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SUBROUTINE    PREDBSCEI , E J,NP I , NPJ, NC P) 
IMPLICIT    REAL*8< A-H, P-l) 
REAL*4    Y(450),YS (50),DM< 50) 
DIMENSION    X<450)  ,XS(50) 

COMMON    TENO, YMAX, THRESH, EPS, X,XS,Y,YS, DM, NL,NR,NK 

DATA    IWR/O/ 
C 

C       KNOT     PPEnTfTtPN     BY    PINAFY    SEARCH    (     BISECTION    ): 
r  NP=(NPI+NPJ)/2  SINCE    NPI     <    NR     <    NPJ. 

C 
10  NP=(NPI*NP,J  )/2 

NF=NR+1 

CALL    PCUBE(NE,E*AX) 

H=X(NR)-X(NL ) 

MCP  =  Nrp-KMR-NL) 

I»=(T  WP.NF.O) 
♦WPTTE(6,10r>  )    NL,NP,  X  (  NR  )  ,  Y  (  NR  )  ,  H,  EMAX 
100       FORMfiT(i0X,2(2X,I3) ,4 ( 2X  ,E12. 5 ) ,3X, 1H*» 

IE(EPS-EMAX)     30,20,40 
?n  EI=PPS 

PETUPN 
30  EJ=EMAX 

NPJ=NR 

00    TP    50 
40  EI  =  FM*x 

NPT=NR 
C       CUCUL*TE    THE    NEXT    KNOT. 
50  IFCNPJ-NPI     .LT.    6)       GO    rp    60 

GO   Tf    10 
r       IF    nHERE    ARE    LESS    THAN    6    PTS.    IN    BETWEEN,    RETURN. 
60  EI=^m*x 

RETURN 

END 
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